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ABSTRACT 


A program designed to numerically integrate the equations of motion of an 
artificial earth satellite is described. The method used is the Cowell predictor- 
corrector process with a local error control which is based on a variation of 
order and/or stepsize during the integration, as described in Velez, C. E., 
"Local Error Control and its Effects on the Optimization of Orbital Integration," 
NASA X-553-67-366. The effects of such controls on the integration are dis- 
played in the form of a resume which includes computer timing data, a stepsize 
distribution analysis, and a propagated numerical error estimate. In addition, 
an option to numerically integrate in multirevolution steps is described. 
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I. DESCRIPTION 


In the computation of orbits of artificial earth satellites by numerical inte- 
gration, advances in the theory of perturbations have led to the use of more 
complex force models and have resulted in the need for more efficient integra- 
tion techniques. On the other hand, the availability of highly accurate tracking 
systems has increased the need for correspondingly accurate numerical 
processes. 

The program to be described is the result of an investigation into numerical 
integration methods as applied to the calculation of orbits. In particular, the 
program is designed to examine how controlling the local error affects the effi- 
ciency and accuracy of the process. 

The numerical integration method considered is the multi-step technique 
using the well-known "summed" form of the Adams-Cowell formulas.* During 
this process the local error is controlled (restricted) by variation of the parame- 
ters p and h, the order and stepsize associated with these formulas. Variation 
of the stepsize and/or order is governed by specifying a set of error tolerances 
or upper and lower bounds on the order. (See Section IV.) 

The effects of such control on the integration are displayed in the form of 
a resume which includes timing estimates, a stepsize distribution analysis and 
a propagated numerical error estimate.** As an application, these results 
could be used in the calibration of other numerical integration techniques which 
require the specification of p or h. 

An additional feature has been incorporated into the program which allows 
one to use a multirevolution process along with the integration, when many revo- 
lutions of orbit are required. This process is described in Section V. 


*EquivaIent Methods: (1) Adams-Moulton; Stormer-Cowell. (2) Gauss-Jackson (central differences). 

**For the case of 2-body (elliptic) motion. 
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II. FORMULAS FOR INTEGRATION MODEL 


The integration model employed is a multi-step process using the following 
form of the Adams-Cowell formulas x : 


(Predictor) X n+1 = h 2 ^ XI S n + ^ X n + ^VX + ^ V 2 X +.... + a k V k ‘ 2 X (1) 


(Corrector) X n + 1 = h 2 [ n S„ +^X n+1 - ^V 2 X n + 1 + .... + a* V k " 2 X n + 1 J 


(1)' 


(Predictor) 


X n+1 = h[^S n+ i-X +1 |vX + .... + a k V k -iX n ] 


( 2 ) 


(Corrector) X n+1 = h pS„ + 4 x „ + i “ ll vx n + i + + < vk_1 x „ + i J 


(2)' 


where 

X = (X, Y, Z), X = (X, Y, Z), h is the step size and p = k + 1 is the order; the 
n S n , I S n are the 2nd and 1st "sums" respectively. 

The coefficients — cr. , o-*, a., a* — are obtained by using the recurrence 
relationships: 


<T 


0 


1 


a 


m 



J 

j + i 


( 3 ) 


iHenrici, P. (1962): "Discrete Variable Methods in Ordinary Differential Equations”, John 
Wiley & Sons, Inc., New York, pp 192-195 and pp 291-293. 
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The backward differences are computed using the relationship 


V ' X = y n-i 


V" 


z__, 

i = 0 


( 4 ) 


The first and second "sums" are defined by 


V " 1 X = X S , V " 2 X = v 


n n 


1 {%) - ll \ 


and are computed by inverting the corrector formulas for X and X : 


ir 


X i .. i 

; J_ y _ 

Vi h 2 [12 » 24 


240 v X n + • • • + CT k 


* yk-2 x 

k n 


*S = -T- ~ 
n 1 h 


[t \ - i 


( 5 ) 


12 VX „+ + < Vt '‘ X 
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then updated using 


*S n = I S. + X 

n n “ 1 n 


n s = u s . + *s . 

n n - 1 n 


( 6 ) 


Remark: In order to keep 1 *S the same order of magnitude as X, the equations 
to be integrated are taken as h 2 X. The above formulas were modified 
to allow this adjustment. 
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HI. COMPUTATION PROCEDURE 


3.1 INTEGRATION STARTING PROCESS 

Formulas (1) and (2) required (k) starting values of the solution in order to 
obtain the (k - 1) backward differences required for the given order (p). 

The starting values (See Figure 1) are computed by employing high order 
Runge-Kutta type formulas (Appendix A) which are particularly suited because 
of the degree of accuracy obtainable. The computation procedure is then: 

(1) Given the initial position* vectors X o = (X o , Y o , Z o ), the order (p), and 
the step size (h), the values X. = X( t + ih), i = 1, 2, .... k - 1 are 
computed. (The velocity vectors X. = (X if Y. , Z.) are also produced 
by this process.) The acceleration vectors X. = (X i9 Y., Z. ), i = 1, 

2, .... k - 1 are computed using the equations of motion given in Ap- 
pendix B. 

(2) Using (4) with r = i , the backward differences V* X n , i = 1, 2, .... k -1 
are produced, and using (5) and (6) with n = k - 1, : S n and XI S n are 
obtained. 

Once the starting table is complete, the integration proceeds as in Sec- 
tion 3.2. 


3.2 COWELL INTEGRATION PROCEDURE 

Having generated a table of starting values, (information above the line in 
Figure 1), the Cowell integration proceeds as follows: 

(a) Compute X£ + 1 (n = k -1), using the predictor formula in (1). 

(b) Compute the backward differences V 1 X n + 1 , i = 1, 2, .... k - 1 as 
follows: 


— — — •• •• 

For simplicity in describing the procedure, only X - (X, Y, Z) and X - (X, Y, Z) will be used. 

* • • • 

Also note that no further references to X ~ (X, Y, Z) will be made since these can be obtained 
by performing the operations described for integrating X, 
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VX n + 1 = X + 1 -X n 

V2 ^n + 1 = V ^n + 1 


V k_1 x +1 = V k_2 x- , - V k ~ 2 X . 

n + 1 n + 1 n 


(c) Obtain a corrected value for X* +1 using (1)' . 

(d) Test |Xp +1 - X^ +1 ! < S , where (S) represents a specified allowable 
tolerance,. Since X = (X, Y, Z), if any one of the three coordinates fails 
the test, X n+1 is evaluated and steps (b), (c), and (d) are repeated with 
Kii replacing XP + 1 and X c ^} replacing X« +1 in (d). 

(e) If the difference of two successive values of the solution in step (d) < S , 
the iterative process is complete and the retained X n+1 is for the former 
value of X n+1 . (Note: If the first corrected value of X +1 is acceptable, 
the step is complete with only one evaluation of the derivative.) 

Since two corrections are usually sufficient 2 , the integration is termi- 
nated if the number of iterations exceeds three. It might then be ad- 
visable to change the order (p) or the step size (h) or if the prescribed 
predictor-corrector tolerance requirements are not too stringent, the 
tolerance (S) might be increased to acquire a more rapid convergence 
of the predictor-corrector cycle. 

(f ) This completes the predictor-corrector cycle for the evaluation of a 
point. The index (n) is updated by one, the sums are updated using (6), 
and the integration process is continued by returning to step (a). 

As the integration proceeds the requested vectors are supplied as output 
data. If a requested vector is not one generated by the integration process, it is 
produced by interpolation. 


2 Hull, T. and Creemer, A. (1963): "Efficiency of Predictor-Corrector Procedures”, J. ACM, 
Vol. 10, pp. 291-301. 
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IV. LOCAL ERROR CONTROL 


An improvement in the efficiency of the integration process is obtained by 
controlling the local truncation error 3 . Associated with the integration formula 
(1) is an approximation to the local error given by 


U = 


h 2 VP -3 X I 

1 k n 1 


(?) 


the magnitude of last difference vector retained in the formulas. 

It is evident from (7) that controls on the local error are directly dependent 
on the variability of the stepsize (h) and the order (p). The program has, there- 
fore, been designed to alter these parameters during the integration as U n varies. 
In particular, the order and/or stepsize can be varied so that the relation 
Tj > |uj > T 2 is satisfied for each n, i.e., at each step of the integration, U n 
is tested. If | U n | < T 2 , the order (p) can be decreased or the stepsize (h) in- 
creased. If | U n | > T 1 , the order (p) can be increased or the step (h) decreased. 

4.1 PROCEDURE FOR CHANGING STEPSIZE 

A stepsize change during the integration requires k- values of the solution 
at the new stepsize to provide the required table (Figure 1). Since, in general, 
many values of the solution at the old stepsize are available, an interpolation 
scheme over these values is used to obtain the corresponding values in incre- 
ments of the new stepsize. We remark that if an increase in stepsize is re- 
quired, the last point computed is accepted, and if the step is to be decreased, 
this point is rejected since the associated local error is larger than the allowable 
upper bound T 1# 

In the program, two techniques designed for changing the stepsize are 
available: 

(a) Halving-Doubling: 

When using this option, the stepsize is modified by either halving, or 
doubling the current stepsize. In the case of doubling, the required back points 


3 Velez, C. E. (1967): ’'Local Error Control and Its Effects on the Optimization of Orbital 
Integration”, NASA X-553-67-366. 
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are obtained by selecting every other point of the last 2k points. In the case of 
halving, the last k/2 points are used along with midpoints obtained by interpola- 
tion. 


(b) Stepsize Computation: 

When using this option, the stepsize is computed as a function of the current 
order and local error estimate, i.e., given an "allowable" local error cr 1 , [cr 1 e 
(T 2 , Tj)] , the new stepsize h opt is given by 


h 


Op t 



l/p + 2 


( 8 ) 


If U n < o-j, the stepsize is increased, and if U n > cr 1 , the stepsize is decreased, 
where the new stepsize is approximately "optimal" with respect to a 1 , i.e., the 
"largest" stepsize which will allow the local error cr 1 for the given order p. 

Once the new stepsize has been computed, the required k values are ob- 
tained by interpolation over the backpoints available at the old stepsize. 

In both cases, once the required backpoints at the new stepsize are obtained, 
the computation proceeds starting with step 2 of Section 3.1. 


4.2 PROCEDURE FOR CHANGING ORDER 

The process of changing the order is less involved than that of changing the 
step. Because of the type of formulas being used, changing the order amounts 
to decreasing or increasing the number of terms retained. If the order is to be 
decreased, one less difference is retained in the computation of subsequent 
points. If an increase in order is required, the last point is rejected, the order 
is increased, and the point is recomputed by returning to step 2 of Section 3.1. 

Since the local error is allowed to vary through a range, the order being 
used is one that will satisfy the tolerances but will not necessarily be the 
smallest or "optimum" order that can be used. An option has been included to 
test the differences at each step against a tolerance cr 2 until the optimal order 
corresponding to this tolerance is established. Although T 2 < a 2 < T x , cr 2 
should be close to T x to obtain a "smallest" order whose local error will satisfy 
the upper bound. 
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4.3 MECHANICS FOR VARYING ORDER AND STEP 


Depending on the nature of the orbit being integrated, the step and order 
can be varied alone, in combination*, or can. remain fixed throughout the inte- 
gration. These operations are referenced as (a) vary order - fixed step, 

(b) vary step - fixed order, (c) vary order - vary step and (d) fixed order - fixed 
step modes. 

The mechanics of the operations are as follows: 

(a) The controls for vary order - fixed step are applied such that for 

U n < T 2 , the order (p) is decreased by one; for U n > T x , the order is 
increased by one. 

(b) The controls for vary step - fixed order are applied such that for 
U n < T 2 , the stepsize is increased; for U n > T lf the stepsize is de- 
creased. 

(c) The controls for vary order - vary step are applied such that for 

U n < T 2 and p < L x , the stepsize is increased; for U n > Tj and p > L 2 , 
the stepsize is decreased, where Lj and L 2 are integers specifying the 
lower and upper bounds on the order, i.e., the stepsize is modified if p 
fails the inequality L 2 >. p > Lj . The order (p) is then arbitrarily ad- 
justed to Lj + (L 2 - Lj)/2. After one point is computed at the new 
stepsize, (p) is optimized. 

If, however, U < T„ and p > L,, the order is decreased; if U > T, 
with p < L 2 , the order is increased. 

The modes (a), (b), and (c) when selected are also used to adjust the initial 
starting table. Mode (a), by varying the order, insures that the best order for 
the given step is selected. When modes (b) or (c) are selected, the table is re- 
started if the stepsize requires changing. 


Both step and order optimization can be used in combination with a restriction cr j > cr 2 . 
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V. THE MULTIREVOLUTION PROCESS 


5.1 FORMULATION 

An option has been included to "step ahead" the calculations in multirevolu- 
tion increments. The multirevolution procedure is a combination extrapolation 
and integration technique. 

A cycle in the algorithm consists of first extrapolating or "predicting" the 
orbital elements n - revolutions ahead. Then, starting with these extrapolated 
values, the equations of motion are integrated over one revolution. Finally, the 
extrapolated values are improved by using a corrector formula along with the 
information obtained from the single revolution integration. 

Consider the following definitions as they relate to the formulas and pro- 
cedure to be outlined: 

(a) f ■ — the value of an orbital element at the descending node of the j th 
revolution. Since this is a 3 -coord mate system, every reference to a 
function actually represents a computation for each coordinate, both 
position and velocity. 

(b) n — the number of revolutions to be stepped ahead. 

(c) V — the backward difference taken at n times the step of V, so that 
ifVf j = fj - f j - 1 * then V f. = f. - f._ n . 

(d) k — the order of the highest backward difference to be retained in 
the multirevolution formulas. 

The formulas used in the multirevolution process 4 * are: 

k 

(Predictor) A n f. = n V I ( Af i) < 9 > 

i = 0 


4 Velez C. (1967): Notes on the Numerical Integration of Orbits in Multirevolution Steps, NASA 

X-542-67-341. 
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(Corrector) 


V f. 

n J 


K 

n L< v -( M .) 


i “ 0 


( 10 ) 


where 


Af. 

j 


f i + 1 - f i 


1 - 


2Z b i “i-j 

j = l 


i = 1,2,3. 


a* = - / b. a* . 

1 / . > 
i = i 


i = 1,2,3.. 


( - 1 > k *‘ ^(4* i) 


i_=l 

(k + "l)“ ! 


k = 1,2,3, 


where 


a 


0* 


a O> t > 0 


1 


5.2 COMPUTATION PROCEDURE FOR MULTIREVOLUTION ALGORITHM 

Just as in the Cowell integration, the multirevolution formulas require a 
set of starting values. This starting table (see Figure 2) is obtained by inte- 
grating using the Cowell procedure until a sufficient number (kn+ 2) of values 
of orbital elements at the descending node of a revolution are obtained. (Values 
of the orbital elements at the descending node of any revolution are obtained by 
inverse interpolation, i.e., a direct interpolation is made to find the time of 
nodal crossing, followed by an inverse interpolation to obtain the values of the 
orbital elements.) Then the forward differences Af. = f. + 1 - L, i = 0, n, 

2n kn are computed, and from these values the differences V‘ (Af kn ), 

i = 1, 2, .... k; where 
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Orbital Elements 1st Backward 

at Descending 1st Forward Diff. at n Times 


Node Difference the Step 



Values above the line are required starting information. 

Those below the line are then produced using the extrapolation procedure. 


Figure 2 
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(k-l)n 


(ID 


V 

n 


( Af kn) 


Af kn * Af 


V n (Af kn ) = Af kn - 2Af (k . 1)n + 


Af 


( k - 2 ) n 


Having obtained the necessary starting information, the extrapolation cycle 
can begin: 

(1) Extrapolate n -revolutions ahead using (9) with j = kn to obtain A n f kn . 

(2) Compute a predicted (extrapolated) value f (k + 1)n , — the element at the 
descending node of the (k + 1) n th revolution by 

f(k+l)n “ Af kn + f( k -l)( n + l) - 

(3) Using the extrapolated values, integrate one revolution to the next node. 
Then compute the forward difference Af (k+1)n = f (k+1) - f (k+1)n * 

(4) Compute the backward differences V n ( f (k + l)n) i = 1. 2 ... k using 
( 11 ). 

(5) Using (10) to obtain a corrected ^ n (f( k +i) n )> a corrected value of the 
orbital elements f (k+1)n = V B (f {k+1)n + f (k - 2)n+1 )ls computed 
and the cycle is complete. 

The multirevolution process is continued by repetition of the 5 cycle steps 
with (j) in step 1 updated for the current revolution, i.e., j = (k + l)n, (k + 2)n, 


Remark: There is an option in the program to delete step (5), i.e., to extrapolate 
using only the predictor formula. 
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VI. FLOW CHARTS 



Flow Chart A 
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ERROR CONTROL 


1 REQUIRED BACKPOINTS FOR A STEP CHANGE ARE COMPUTED USING INTERPOLATION. 


Flow Chart B 
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Flow Chart C 
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VII. OPERATING INSTRUCTIONS 


7.1 SEQUENCE OF INPUT DATA CARDS 

1. Logical switch settings 

2. Cowell Integration beginning and ending times 

3. Step sizes for Runge Kutta integration 

4. Epoch date 

5.-6. Orbital elements at epoch 

7. Tolerances for controlling local error 

8. Mode of Operation 

9. Cowell integration order, step or order limits. 

10. Special print out interval 

11. Extrapolation mode — order and stepped revolutions. 

During a computer application in which several integrations are to be per- 
formed, data cards 1, through 9 are required for the first integration while suc- 
cessive integrations require only the logical switch settings (Card 1) and changed 
information. 


7.2 DESCRIPTION OF INPUT CARDS 


Card #1 — Format (30L1, II) A true (T) in the column corresponding to the 
ISWT index number indicates: 


Column Format 


Remarks 


1 


T ISWT(l) — New run — Read all new parameters. 


2 T ISWT (2) — Repeat last run — Different mode. 

3 T ISWT (3) - Repeat last run — Different orbital elements. 
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Column 

Format 

Remarks 

4 

T 

ISWT(4) — Repeat last run — Different tolerances. 

5 

T 

ISWT(5) — Repeat last run — Different order or step. 

6 

T 

ISWT(6) — Repeat last run — Different integration times. 

7 

T 

ISWT(7) — Use with options other than 2 or 5 to suppress 
print out except for a specified portion of the orbit. 

8 

T 

ISWT(8) — Use NPT to produce print out based on the 
number of integration steps. 

9 

T 

ISWT(9) — Used to indicate that the initial conditions 
are to be input as position and velocity vectors. 

10 

T 

ISWT(10) — Repeat last run — Change Runge-Kutta step- 
sizes. 

11 

T 

ISWT(ll) — Run in extrapolation mode. 

12 

T 

ISWT(12) — Print error vectors for elliptic motion. 

13 

T 

ISWT(13) — Restart with new extrapolation case. 

14 

F 

ISWT(14) - Not used. 

15 

T 

ISWT(15) — Extrapolate by prediction only. (False — 
by prediction and correction.) 

16 

T 

ISWT(16) — Run with step optimization. 

17 

T 

ISWT(17) — Run with order optimization. 

18 

T 

ISWT(18) — Print out nodal information. 

19 

F 

ISWT(19) — Not used. 

20 

F 

ISWT (20) — Not used. 

21 

F 

Full earth gravity applied. (True -elliptic motion only.) 


19 


Column Format 


Remarks 


22 

23 

24 

25 
31 


T Lunar gravity applied. 

T Solar gravity applied. 

T Atmospheric drag applied. 

T Solar radiation applied. 

I IND — Integer ^ 0, program execution will be terminated. 


Card #2 — Format (2D15.8,I8) 


Column Format 


Remarks 


1-15 

± X.XXXXXXXXD ± XX 

16-30 

± X.XXXXXXXXD ±xx 

31-38 

mum 


T 0 — Integration starting time (minutes). 

FT — Integration final time (minutes). 

NPT — Print out interval. For normal 
output [ISWT(8) = false] NPT is inter- 
preted as minutes of orbit. If [ISWT(8) 

= true] , NPT is interpreted as the num- 
ber of integrated points. 



Card #3 

— Format (2D10.3) 

Column 

Format 

Remarks 

1-10 

±X.XXXD±XX 

HI — Runge Kutta step size when inte- 
grating forward (c.u.t.). 

11-20 

±X.XXXD±XX 

H2 — Runge Kutta step size for inte- 
grating backward. 




Card #4 — Format (I6,I4,D7.4) 

Column 

Format 

Remarks 

1-6 

YYMMDD 

Year, month, day of epoch. 

7-10 

HHMM 

Hours and minutes of epoch date. 

11-17 

xx.xxxx 

Seconds of epoch date. 
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Cards #5-6 — Format (3D24.17/3D24.17) 


Column 

1-24 

25-48 

49-72 

1-24 

25-48 

49-72 

1-24 

25-48 

49-72 

1-24 

25-48 

49-72 

Column 

1-10 


Format 

± X.XXXXXXXXXXXXXXXXXD ± XX 

± X.XXXXXXXXXXXXXXXXXD ± XX 
± X.XXXXXXXXXXXXXXXXXD ± XX 
± X.XXXXXXXXXXXXXXXXXD ± XX 

± X.XXXXXXXXXXXXXXXXXD ± XX 

± X.XXXXXXXXXXXXXXXXXD ± XX 
— or — 

± X.XXXXXXXXXXXXXXXXXD ± XX 

± X.XXXXXXXXXXXXXXXXXD ± XX 
± X.XXXXXXXXXXXXXXXXXD ± XX 
± X.XXXXXXXXXXXXXXXXXD ± XX 

± X.XXXXXXXXXXXXXXXXXD ± XX 
± X.XXXXXXXXXXXXXXXXXD ± XX 


Remarks 


ELEM(l) — Semi-major 
axis (c.u.l.). 

ELEM(2) — Eccentricity. 

ELEM(3) — Inclination (rad.). 

ELEM(4) — Mean anomaly 
(rad.). 

ELEM(5) — Argument of 
perigee (rad.). 

ELEM(6) — Longitude of 
node (rad.). 

Xl(3) — X 0 , Y 0 , Z 0 , coor- 
dinates for initial position 
vector, (c.u.l.) 


XD1(3) - X 0 , Y 0 , Z 0 coor- 
dinates for initial velocity 
vector, (c.u.i./'c.u.t.) 


Card #7 - Format (5D10.3) 


Format 


Remarks 


±X.XXXD±XX TOL1 — Upper tolerance limit for local trunca- 
tion error control. (T 1 ) 
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Column 

Format 

Remarks 

11-20 

± X.XXXD ± XX 

TOL2 — Lower tolerance limit for local trunca- 
tion error control. (T 2 ) 

21-30 

± X.XXXD ± XX 

CTOL — Accuracy requirement for integration 
corrector cycle convergence. (S) 

31-40 

±X.XXXD±XX 

TOL3 — Tolerance used for optimum step com- 
putation. (CTj) 

41-50 

± X.XXXD ± XX 

TOL4 — Tolerance used for optimum order 
option. (ct 2 ) 



Card #8 — Format (11) 

Column 

Format 

Remarks 

1 

I Mode — 1 - indicates vary step vary order option. 



2 - indicates vary step fixed order option. 

3 - indicates vary order fixed step option. 

4 - indicates fixed order fixed step option. 


Card #9 

- (With Mode = 1) Format (213) 

Column 

Format 

Remarks 


1-3 HI Ll — Lower order limit for vary order-vary step mode. 

4-6 III L2 — Upper order limit for vary order-vary step mode. 


Card #9 - (With Mode = 2) Format (D24.17.I2) 

Column Format Remarks 

1-24 ± X.XXXXXXXXXXXXXXXXXD ± XX DUM — Dummy variable. 

25-26 H ORDER - Cowell integra- 

tion order 
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Card #9 — (With Mode = 3) Format (024.17,12) 


Column 

1-24 

25-26 

Column 

1-24 

25-26 

Column 
1-10 
11-20 
21 - 30 

Column 

1-3 

4-6 


Format Remarks 

± X.XXXXXXXXXXXXXXXXXD ± XX DEL — Cowell integration 

stepsize in minutes. 

II DUM — Dummy variable . 


Card #9 — (With Mode = 4) Format (D24. 17,12) 

Format Remarks 

i X.XXXXXXXXXXXXXXXXXD ) XX DEL — Integration stepsize 

(mins.). 

n ORDER — Order to be used 

for integration. 


Card #10 - (With ISWT(7) = TRUE) Format (3D10.3) 


Format 
± X.XXXD ± XX 
± X.XXXD ± XX 
± X.XXXD ± XX 


Remarks 

BEGTIM — Starting time for print out (mins.). 

ENDTIM — End time for print out (mins.). 

BEGT2 — Start print out at this time and con- 
tinue to end of run (mins.). 


Card #11 - (With ISWT(ll) = TRUE) Format (213) 

Format Remarks 

HI NEXT — Number of revolutions to be ’’stepped". 

HI KEXT — Number of differences retained in the extra- 

polation formulas. 
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7.3 BLOCK DATA 


Input parameters are also supplied through a block data section of the pro- 
gram. These are parameters required by the force model and may be changed 
by replacing the values in the following data statements: 

(1) DATA GM/1.D0/, AE/l.DO/ 

GM — gravitational constant times mass of the earth. 

AE — semi-major axis of reference ellipsoid in c.u.l. 

(2) DATA EMASS (1)/.012299896D0/ 

EMASS (1) — ratio of mass of moon to mass of earth 

(3) DATA EMASS (2)/332951.3D0/ 

EMASS (2) — ratio of mass of sun to mass of earth. 

(4) DATA WE, F, CAPR, CD/.5883369D-01, 298.25D0, 6378.166D0, 2.3D0/ 

WE — angular rotation of earth 
F — earth flattening constant 
CAPR — 1 c.u.l. in km. 

CD — drag coefficient 

(5) DATA AREA, SATMAS/.XXXX D±XX, .XXXX D±XX/ 

AREA — area of satellite in grams/cm 2 . 

SATMAS — mass of satellite in grams 

(6) DATA RHOl, RH02, RH04/2*.XXXD ±XX, 30.D0/ 

RHOl — n 

> differential correction parameters 
RH02 - 1 J 

RHOl — atmospheric bulge angle in degrees 
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(7) DATA ((CS (I, J), I = 1, 20), J = 1, 23)/ 


CS — harmonic coefficients 


25 


Vm. SUBROUTINES USED 


The Cowell integration program is designed as-a system of subroutines 
under the control of an executive routine. The FORTRAN name and brief de- 
scription follow : 


EXEC 

PCCOFF 


TRANS 


TRANS2 


— Directs the integration process. 

— Computes the coefficients for the Cowell integration formulas. 

Calling Sequence: CALL PCCOFF (PX, PXD, CX, CXD, I). 

where I = the number of coefficients to be 
computed. 

— Converts orbital elements to position and velocity coordinates. 
Calling Sequence: CALL TRANS. 

— Converts position and velocity coordinates to orbital elements. 
Calling Sequence: CALL TRANS2. 


TABLE 


FRCS 


CSTEP 


SUMS 


CKDIFF 


— Using Runge Kutta integration, computes the starting table of 
values required for the Cowell integration model. 

Calling Sequence: CALL TABLE. 

— Computes accelerations using the equations of motion (Ap- 
pendix B) or elliptic motion. 

Calling Sequence: CALL FRCS. 

— Performs Cowell predictor-corrector cycle. 

Calling Sequence: CALL CSTEP. 

— Computes X S and n S for Cowell integration. 

Calling Sequence: CALL SUMS. 

— Computes the i th backward differences V 1 X n for Cowell 
integration. 

Calling Sequence: CALL CKDIFF (I). 
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TEST 
TABLE B 
HE MINT 
RK 

OUTPUT 

RESUME 

TNODE 

EXTCF 

FKDIFF 

EPHEM 

EPHQAN 

EGRAV 

SLGRAV 

DRAG 


— Test and control section for local truncation error associated 
with Cowell integration. 

Calling Sequence: CALL TEST. 

— Computes the new starting points when a change of step size 
is required. 

Calling Sequence: CALL TABLEB. 

— Hermite interpolation used to produce backpoints when changing 
step size or to produce points at print out request times. 

Calling Sequence: CALL HEMINT. 

— Performs Runge Kutta type integration using formulas in 
Appendix A. 

— Formats printed output at print request times. 

— Collects information throughout the integration to produce an 
efficiency summary. 

— Extrapolates a point n-re volutions ahead. 

— Computes coefficients for extrapolation formulas. 

— Computes the differences V* of the larger forward differences 
for extrapolation formulas. 

— Computes lunar and solar ephemerides. 

Calling Sequence: CALL EPHEM. 

— Obtains the lunar and solar quantities for a 5-day interval and 
obtains the coefficients for a least squares fit to a 4th order 
polynomial . 

— Computes the acceleration due to earth's gravity. 

— Computes the acceleration effect due to solar and lunar gravity. 

— Computes the acceleration effect due to atmospheric drag. 
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IX. CODING INFORMATION 


BEGT2 

BEGTIM 

CAPR 

CDEL 

CDP 

/COFFS/ 

COWL 

/COWS/ 

CSAVE 

C 1 
cc J 

CTOL 

cx ^ 
CXD J 

DEL 


Defined Symbols* 

— start printout at this time and continue printing to end of run. 

— begin printout at this time and print to a specified final time 
(ENDTIM). 

— conversion factor meters/c. u.l.** 

— Cowell integration step size h in c.u.t.** 

— atmospheric drag constant 

— name of common block containing CSAVE. 

— running time for routine CSTEP. 

— common block containing variables SI, S2, PX, PXD, CX, CXD, 
CTOL, SAVE, ITER. 

— location for saving position predictor and corrector coefficients 
for last retained term of series. 

— predictor -corrector coefficients for extrapolation process. 

— accuracy criterion for Cowell corrector cycle convergence. 

— arrays of coefficients for Cowell corrector formulas. 

— Cowell integration step size in minutes. 


Symbols relating to the constants used in the force model have been omitted; see BLOCK DATA 
subroutine. 

* ^Canonical unit of length (c.u.l.) = 6378.166 kms. 

Canonical unit of time (c.u.t.) - 806.81242 secs. 
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DIFF 

ELEM 

/EMS/ 

ENDTIM 

FNP 

FORCS 

FT 

FX -\ 
FXD > 
FXDD 

hi 

IE NT 

INCOWL 

/INTERP/ 

ISCT 

ISWT 

ITER 

KEXT 

K1 


— array containing differences V 1 X n to be used in the Cowell 
equations for correcting. 

— array containing orbital elements. 

— common block containing E LE M. 

— print out end time to be used with (BEGTIM). 

— used with print request to adjust interpolation location. (FNP = 
2 will cause interpolation between 2nd and 3rd end points of data 
array.) 

— running time for routine FRCS. 

— the total length of integration in mins. 


— arrays for position, velocity and acceleration coordinates inter- 
polated using the Hermite interpolation subroutine. 


— Runge Kutta step size for forward integration. 

— number of step sizes selected during the integration. 

— number of integration steps. 

— common block containing FX, FXD, FXDD, Tl, Kl, M, Ml 
(Parameter Ml in HERMITE routine is called L.) 

— counter for the number of points at a particular step size. 

— set of logical switches to control program operation. 

— number of Cowell corrector cycle iterations. 

— the order to be used for extrapolation. 

— starting point of array to be interpolated using Hermite inter- 
polation. 
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K 


LI 

L2 


} 


Ml 

MODE 


M 

NCT 

NL1 1 
NL2 J 

NEXT 

/NODE/ 

NPT 

N 

/ODEL/ 

/OPT/ 

ORDER 

PERIOD 


location of current integrated point(s) in the vector 
arrays. 

limits on order to be used with vary order — vary step mode 
to control step size changes. 


/LIMITS/ — common block containing TOL1, TOL2, TC, ISCT, ISWT, SW, 
ORDER, LI, L2, MODE. 


— location of interpolated value returned from HE MINT routine. 

— integer value to indicate the manner in which the program is 
operating. 

1 = Vary Step — Vary Order 

2 = Vary Step — Fixed Order 

3 = Vary Order — Fixed Step 

4 = Fixed Order — Fixed Step 

— order of hermite interpolation polynomial. 

— total number of Cowell steps taken. 

— maximum and minimum order limits. 

— number of revolutions to be "stepped" in extrapolation process. 

— common section containing XNODE, XDNODE, TNOD, C, CC, 
KK, Jl, NDIFF, NEXT, KEXT, INODE. 

— output interval. 

— number of differences carried in the integration. 

— common block containing NL1, NL2. 

— common block containing BEGTIM, ENDTIM, BEGT2. 

— order of the Cowell integration formulas. 

— computed period of the satellite. 
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/PERTS/ 

PXD 
PX 
SAVE 


51 

52 


} 


sw 

T 

TC 

TEPOCH 

/TIMING/ 

TI 


TNOD 

TT 

TREQ 

TOL1 

TOL2 

TOL3 

TOL4 

TOUT 

TSEC 


— common block containing ALT, DEN, HM, WE, ESQ, B, F, CDP, 
RHOl, RH02, RH04, CAPR, AREA, SATMAS, CD, MD. 

— arrays of coefficients for Cowell predictor formulas. 

— array containing differences V 1 X n to be used in the Cowell 
equations for predicting. 

— arrays containing first and second sums, I S, n S. 

— set of logical switches for internal program control. 

— elapsed time from epoch in c.u.t. 

— constant for converting minutes to c.u.t. 

— integration starting time in minutes. 

— common block containing stored information for a print out 
resume. 

— interpolation time for positions and velocities when using 
Hermite interpolation. 

— array containing times of nodal crossings. 

— square of Cowell integration step size in c.u.t. 

— output request time in minutes from epoch. 

— upper and lower bounds on local truncation error. 

— crj , ct 2 defining "allowable" local error for step and order 
computation. 

— elapsed time from epoch in minutes. 

— conversion factor sec/cut. 
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/WORKER/ 

XDD 

XDNODE 

XD 

XNODE 

X 

XXDD 


— common block containing X, XD, XDD, XXDD, DIFF, CDEL, TT, 
T, TREQ, TOUT, K, N. 

— array for storing acceleration coordinates of the form h 2 X. 

— array for storing extrapolated velocities. 

— storage array for integrated velocity coordinates. 

— storage array under extrapolation mode for position coordi- 
nates. 

— storage array for integrated position coordinates. 

— value for current acceleration coordinate computed without 
h 2 in the equation. 
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X. ERROR CODES 


INITIAL TABLE FAILED TOLERANCES PROCEDURE 

RESTART WITH NEW STEPSIZE. 

RUN TERMINATED . . . ITERATIONS GREATER THAN 3. 

FINAL TIME = .XXX D±XX 

ORDER CYCLE DOES NOT CONVERGE WITHIN SET LIMITS. 
FINAL TIME = .XX D±XX 

STEP CHANGES GREATER THAN 90. RESUME INCOMPLETE. 
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XI. STORAGE REQUIREMENTS (1108 configuration) 


INSTRUCTION SPACE REQUIRED 


No. of Decimal Locations Approximate 


EXEC 

CKDIFF 

COEFF 

CSTEP 

DIFF 

DNVERT 

DRAG 

EGRAV 

EPHEM 

EPHQAN 

EXTCF 

FIT 

FKDIFF 

FRCS 

HE MINT 

OUTPUT 

PC C OFF 

RESUME 

RK 

RYMDI 

SLGRAV 

SOLRAD 

SUMS 

TABLE 

TABLE B 

TEST 

TNODE 

TRANS 

TRANS2 

YMDAY 

NTAB$ 

SHADOW 
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11.2 


DATA SPACE REQUIRED 


Labeled Common 


No. of Decimal Locations Approximate 


3978 

61 

14 

4 

879 

2886 

556 

1085 

4 

2 

2 

4 

4 

6 

3 

20 

8 

29 

14 

11 

96 

13 

922 

433 

96 

6 


Additional Storage: 


Approximate Locations 


Parameters and constants 
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XII. PROGRAM LISTING 
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ooooooooe© 


£LT EXEC.1.671009. 48 S38 

EuF in 


COWELL ORBIT ERROR AND TIME ANALYSIS PROGRAM 


DOUBLE PRECISION X » XJ. XOD.DIFF , P A »PXD»CX . L XU » T . TC » DEL . COEl » TO- iT . Si 
1 » S2 . XI . Xul . XXUD. TEPOCH > TOL1 » TOL2 » CTOL » FT » TT . DSQRT . HI » H2 » El EM. 
2 ENDTIM » BEGT I i*i » FT l » PERIOD » CSAVE . STEPO » O loT . SAVE » HUM » TM » SDEl . S 
30RD » TEMP3 . TEMP4 . 8tGT2 , 7 0L3 . T0L4 
DOUBLE PRECISION F X. FXD.FXDO » Tl » TREQ . F.JPT .FNP . RSAVE » SFNPT , TER 
DIMENSION XI (3) . XDl (3) 

DOUBLE PRECISION xNODE.XDNOOE.TNOD.C.CC 

DOUBLE PRECISION 7 HETGo . 7HDOT1 , Thl> 0T2 .OAYl » QMlN.ETIME , 

* ORAO»i.T*OPI .O kSEC.GM. AE.R.RSO.RQ.THETG. 

* CS.CEi TER.EPSEC.OSTART. YNDAY.TC1.DL0G 

double PRECISION ALT » ULN » HM » ViE . ESC » 8 * RHOl . RH02 . RH03 * RH04 » CN » C P* 

* F.CMPI--,AREA»SATv;AS.CD 

DOUBLE PRECISION xYc.E ASS. PSUn.CSUBR. sigma. TSEC 
integer ORuER 
LOGICAL ISwT.SW 
COMMON/ EMS/ELEM ( 6 ) »7EP 

C0Mm 0N/W0RKEK/X(20U.3> ,XD(200.S) ,xUO( 200.3) » XXDO ( 3 ) » 0 IFF ( NO . 3 ) . 

1 coel.it.t.treq.tout.k.n 

COMMON/LIMI TS/T0L1.T0L2.TC. I SCT. I SwT( 30) »SW(20) .ORDER. LI. I 2. MODE 
COMMON/OPT I M/T0L3.T0L4 

COMMON/ COWS/SI ( 3 ) »Sc(3> »PX(b3 ) .Pxd( 63) .CX(63) »CXD(63) .CTOl »SAVE(60 
1 » 3) » I7El< 

COMMON/NODE/XNOUE (200.3) ,XUt':OUE(2U0.3) .TNOQ(200) . C ( 20 ) . CC ( 20 ) . 
iKK.Jl.NOIFF. NEXT » kExT . INDUE 
COMMON/TIMING/STEPD (90 , 3) . TFMP3* TEMP4 

* .0LUT.TM.XM.YM.2M.XS. I ENT . I7EKS. INCOWL. ICH 
COMMON/INTERP/FXOSO.3) . FXl) (60.3) , F XUU ( 60 » 3) » T I » Kl . M » Ml 
COMMON/COFFS/CSAVt ( E ) 

COMMON/OUEL/NLl . NL2 
COMmON/RKST /N l . H2 
COMMON/OPT/BEGTIM.ENUT IM.BEGT2 
COMMON/ TIMES/TAB.COwL.FORCS 

common/brief/rsave ( 10 ) 

COMMON/CONST/DR AD . QT rfOPl » DRSEC . R AU . RESEC 
COMMON/CONST 1/THETGu ( 111 ) , THDOT 1 . THD0T2 , DMIN. ETImE . 1 YHEG 
C0MM0N/C0NST2/GM. AE.K.RSQ.RD.ThETG. TCI 
C0MMON/CONST3/EMASS(2) ,XY2<4) .LS 
C0MM0N/C0NST4/PSUN . CSUhR . SIGMA . F 

COMMON/PERTS/ ALT ( 41 ) . OF N ( 41 . 2 ) . HM ( 40 » 2 ) . WE . ESG » 3 . CN » COP. RnOl . kHO 

* RH03. Rim04. C APR. AREA. SATMAS.CD.MD 
COMMON/COFIT/CS (5.9) .QaYI .CENTER. DSTART 
COMMON/CSUN/ ASUN » PMOON , EmOON » T ASUN (10) 

ifa FORMAT (1H0 23X . 2 OHPR 1 GRAM OPERATING IN. 12. 40H REVOLUT 7 OnS EXTR 
♦APOLATION MODE OF ORDER. 12) 

193 FORMAT (213) 

194 FORMAT (3U24. 17) 

195 FORMAT (2D10.3) 

196 F0RMAT(D24. 17.12) 

197 FORMAT ( 5010 • 3) 
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198 FORMAT (2D15. 8, 18) 

199 FORMAT(Il) 

200 FORMAT (30L1, II) 

201 F0RMAT(3D24. 17/3024. 17) 

202 FORMAT (1H1,4BX28HC0WELL NUMERICAL INTEGRATION) 

203 FORMAT ( 1H0 , 49X10HTOLERANCES/1H020X4HTOL129X » 4HT0L229X4HCTnL/1 7X* 
*D10.3,2(23X,D10.3) ) 

204 FORMAT ( 1H040X32HExTRAPoLATION BY PREDICTION ONLY) 

205 • FORMAT (1H033X42HEXTRAPOLATION BY PREDICTION AND CORRECTION) 

211 FORMAT (1H039X31HRUNGE-KUTTA STEPSIZES IN V.U.T./1H040X,3Hh1=D10. 
13X»3HH2=D10.3) 

212 FORMAT ( 1H0 ,40X1 1HF INAL TIME=Dl9.l2) 

213 FORMAT <1HO,45X,7HPERIOD=D19. 12) 

217 FORMAT ( 1H0/ / , 49X * 12H0RDER LIMITS/ 50X. I3,4X*I3) 

244 FORMAT ( 1H *///,18X,40HSTEP AND OROER OPTIMIZATION WITH TOl 3 OF 
* D15.8, 12H AND T0L4 OF D15.8) 

245 FORMAT (1X,///27X,31H0RDER OPTIMIZATION WITH TOL4 OF 015. ft) 

246 FORMAT ( 1H ,///,28X,30HSTEP OPTIMIZATION WITH T0L3 OF D15.ft) 

247 FORMAT (1H1//////////// 26X, 57HVARI-0R0ER VARI-STEP COWELl NUMER 
*AL INTEGRATION PROGRAM) 

248 FORMAT! 1H1/////////////////////////26X58HVARI-0RDER FIXEO-STFP COW 
lELL NUMERICAL INTEGRATION PROGRAM) 

249 FORMAT (1H1/////////////////////////26X59HFIXED-0RDER FIXEn-STEP 
1WELL NUMERICAL INTEGRATION PROGRAM) 

250 FORMAT { 1H1 //// ///// /// //////26X58HFIXED-0RDER VARI-STFP C0U> 

1ELL NUMERICAL INTEGRATION PROGRAM) 

251 FORMAT (1H043X24HINIIIAL ORBITAL ELEMENTS/ 1H016X9HS.M. AXIS24X12HEC. 
1CENTRICITY24X11HINCLINATION) 

252 FORMAT (3 ( 10X ,D24. 16) ) 

253 FORMAT ( 1H015X12HMEAN ANOMALY , 20X15HARG. OF PERIGEE20X13HLONG. OF H 
10DE) 

350 FORMAT (21H0 TIME IN ROUTINE, IX, A6.F12. 5. 4H SEC) 

355 FORMAT (21H TIME IN ROUTINE 1X,A6,F12.5) 

360 FORMAT ( 1H1////////////////////60X11HFORCE MODEL) 

362 FORMAT (///50X,36HA» GEoPOTENTlAL COEFFICIENTS APPLIED) 

364 FORMAT (///50X»18HA. ELLIPTIC MOTION) 

366 FORMAT (///50X»24HB. LUNAR GRAVITY APPLIED) 

368 FORMAT (///50X,24HC* SOLAR GRAVITY APPLIED) 

369 FORMAT (///50X»20HU« ATM. DRAG APPLIED) 

372 FORMAT (///50X,26HE. SOLAR RADIATION APPLIED) 

370 FORMAT (///50X HHEPOCH DATE » 16* IX, 14, IX, 1)15. 4) 

TSEC=806. 8124200 

TC1=TC/8.64D4 
TC=60.D0/TSEC 
C DRAG CONSTANTS 

DO 24 J=l,2 
DO 24 1=1,40 

24 HM(I,J)=(ALT(I)-AlT(I+i) ) / (DLOG (i)EN ( 1 + 1 » J) /GEN ( I , J) ) ) 
CDP=(2.3D0*AREA)/{2.D0*SATMAS) 

F=1.D0/CN 

£SQ=F*(2.D0-F) 

8=(1.D0-F) 

RH04=RH04*DRAD 

C SOLAR RADIATION CONSTANT 

SIGMA= ( CSUBR*PSUN*AkEA ) /SATMAS 
C CONVERT TO CUL/CUT**2 

SIGMA=(SIGMA*TSEC**2)/<CAPR*l.U+05) 

C MAXIMUM ORDER LIMITS 

NL1=4 
NL2=30 
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i\IO=NL2 + l 

THDOT1=THDOT1*ORAU 
00 23 I = lflO 
T ASUN ( I ) =T ASUi J ( I ) *KAL) 

25 THLTGO ( 1 ) =THhTG0 ( 1 ) *Ufi AD 
THD0T2=THD0Tl+DT.v0Pi 
C CONVERSION FROM METEkS TO VOL. 
ASUN=ASUN/(CaPR* 1.0+03) 
PMOON=PMOON/ ( C APR* 1 . 0+i, ;-s ) 

O 0 1 I=1»N0 

CALL PC C OFF (PX»PXG*cX»cXU»I) 

1 CONTINUE 

Sbl KlAU ( 2 » 200 ) I i>«T » IhNu 

IFUENO.NL. 0) GO TO 103 
CALL CLOCK ( XS ) 

TAU=0.0 

C0WL=0.0 

FORCS=0.0 

INCOWL=0 

ITEKS=0 

1CFU0 

IPT = 0 

TLMP3=0.U0 
00 2 I=1»2U 

2 Sh ( I ) =• FALSE • 

00 0 0-1 t 3 

uO 3 I — 1 * 20 

3 STEPO(I» J)=0.u0 
ILNT=1 

C FImP=1NTERP0LA1 ION I OCATION 
FNP=2.D0 

iF(iswna)) fnp=o.du 

IFUSWT(l) ) GO TO 401 
IF USWTC2) ) GO TO 4o4 
IF ( ISWT (3) ) GO TO 4U2 
IF ( ISWT (4) ) GO TO 403 
IF (ISWT (5) ) GO TO 403 
IF (ISWT (6) ) t>0 TO 40l 
IF ( ISWT ( 10 ) ) GO To 423 
IF (ISWT ( 13) ) GO To 425 
4:)1 RLAU(2. 198)TLP0CH,FT»NFT 
FNPT=NPT 
SFNPT=FNPT 
FT1=FT*TC 
TEP=TEPOCH*TC 
IF (ISWT (6) ) GO TO 4u7 
423 RLAU T 2 1 195) H1.FI2 

IF ( ISWT (10) ) GO TO 407 
402 T=TLPOCH*TC 
TOUT=TEPOCH 

RLAU ( 2 ► 192) IEPYrtUt ILPH,v.»EPSEC 
172 FORMAT ( lb » 14 » t)7. 4 ) 

USTART=YMDAY ( lEPYhO * IEPhM • LPSLC ) 

IYBEG=IEPYMD/10000 

IF C ISWT (9) ) GO TO 4u9 

RLAU (2*201) ElEK 

CALL TRANS(X1*XU1) 

GO TO 410 

409 KLAU (2> 201 ) XI * XD1 
CALL TRANS2 ( XI » XU1 ) 
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410 DO 412 1 = 1 * 3 
X(1*I>=X1(I) 

412 XD(1*I)=XD1(I) 

PERI0D=DSQRT(ELEM(l)**3)*6.283185307l795864D0/TC 
IF ( ISWT (3) ) GO TO 411 
4 n3 KEAD (2*197) TOH»ToL2.CTOL*TOL3,TOL4 
IF ( ISWT (4) ) GO TO 407 

404 READ (2*199) MODE 

405 GO TO (413*415*414*416) *MODF. 

413 READ(2* 193) L1*L2 
0RDlR=U1+(L2-L1)/2 
CDEL= 2.D0**(-5) 

N0=2*NL2 

GO TO 417 

414 0RDER=9 

READ (2* 196) DEL 

CDEL=DEl*TC 

N0=NL2 

L1=NL1 

i_2=NL2 

GO TO 417 

415 CDEL=2.0D0**(-5) 

READ (2*196) DUM*ORr>ER 
H=ORDtR 

l2=L1 

N0=2*0RDER 
GO TO 417 

416 KEAD (2* 196) DEL* ORDER 

sw ( b)=. false. 

SW (9)=.FALSE. 

CDEL=QEu*TC 

mo=order 

417 SDEL=CDEL 

sord=ordc;r 

IF ( ISWT (2) . OR. 15W T (b) ) GO TO 407 
IF ( .NOT .ISWT (7) )GC TO 450 
KEAD (2* 194 ) REGTIM* ENDTIM* 8EGT2 
GO TO 450 

425 ISWT ( 11 ) = .TRUE • 

4:j7 DO 408 1 = 1*3 
X(1*I)=X1(I> 

408 XD ( 1 * I ) =XD1 ( 1 ) 

T=TtPOCH*TC 

411 CDEL=SDEL 
ORDER=SOKD 
FNPT=SFNPT 

450 TT=CDEL*CDEL 
D=0RDER-2 
TREG=TEPOCH+FnPT 
IF (ISWT (7) ) TRE0=TREQ+6EGTIFi 
ITT = 1 

IF(FT.GT.TEPUOH) GdTi, (418 *420 , 419. 421 )* MODE 
1TT=2 

fnpt=-fnpt 

IREU=TEPOCH+F -IPT 

odel=-coel 

GO TO (418*420.419*421) .POOL 
416 wK ITE ( 3 * 247 ) 

IF ( ISWT ( 16) .ANu. ISwT ( 17) ) GO TO 44 
IF ( ISWT ( 17 ) ) GO TO 45 
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IF (1SWTU6) ) 
GO TO 40 
/+4 wRITE(3»244) 
GO TO 40 
*5 wR1TE(3.245) 
40 wRITE(3»2l7) 


WR I T E ( 3. ,.'4o) T OL3 

r0i_3» TgL4 

T0L4 
Ll » L2 


4i9 

4 ->0 

4 2l 
422 


WRlTE(3»203)TOLl.T0L2.(.T0L 
GO TO 422 
WRITE13.248) 

I F ( ISWT (17) ) 
wRITE(3.203) 

GO TO 422 
WR1TE(3.250) 

IF' ( 1SWT (16) ) 
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30 


32 

34 

36 

38 

39 

72 

3l 


wRlTt: (3*245) T0U4 
T0L1.T0L2.CT0L 


TOL3 


.. (..RITE ( 3 » 24b ) 

wRITE(3.203) TOL1.T0 l2.cTOL 

GO TO 422 

WK1TET3.249) 

wKlTE(3.203) r0Ll.T0L2.CT0L 
wKlTE(3»251) 

WRITE (3.252) (lLErH I ) » 1 = 1 . 3) 

WRITE (3.253) 

WRITE (3 .252) (lLEM(I) » 1=4.6) 

WRITE(3.213)PLR10o 

WRITE (3.212JFT 

WRITE(3.211) rTl . H2 

IF ( .NOT .ISWT ( 11 ) ) GO To 30 

REA0(2.193)NLAT»Kt.XT 

N0lFF=NEXT*KEXT+2 

DO 30 1=1.3 

DO 35 J=1»I,mOjE 

TNOU(J)=O.OCiO 

XNOl>E(J.I)=0,gO 

XDNODEt J. I )=U.DO 

1N0OE-0 

KK-1 

J1=KEXT+1 
CALL EXTCF 

WRlTE(3.16)t4tXT.KFXT 

IF ( ISWT ( 15) ) wRlTf (3.204) 

IF ( .NOT . ISWT ( lb) ) WRITi- (3.205) 

CONTINUE 

WRITE(3.360) 

IF ( ISWT (21 ) ) GO 10 32 
WRlTE(3.3b2) 

GO TO 34 
WRITE(3.364) 

IF ( .NOT .ISWT (22) ) o 0 TO 36 
WRITE(3»366) 

IF (.NOT. ISWT (23) ) GO TO 36 
«RITE(3»368) 

IF ( .NOT • ISwT (24) ) GO TO 39 
WR1TE(3.369) 

IF (.NOT. ISWT (25) ) GO TO 72 
wRITE(3.372) 

wkITE(3.370) lEPYiol). IEPHR.EPSEC 
K=1 

RSA\/E(4)=T 
OLDT=T 
OAY1=OSTART 
CALL EPHuAN 
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CALL FRCS 
CALL CLOCK (XM) 

CALL TABLE 
CSAVE ( 1 ) =PX ( N+3) 

CSAVE(2)=CX(N+3> 

PX(N+3)=0.D0 

CX(N+3)=0.D0 

RSAVE ( 1 > =DBLE ( FLOAT ( ORDER ) ) 

RSAVE(2>=CDEL/TC 
STEPO (1*1) =CL)EL 
TEMP4=CDEL 
M=ORDER 
TOUT=T/TC 

IF(ISWT(8) ) GO TO 17 
15 SO TO (19»18) » ITT 

19 IF(TREQ.GT.TOUT) SO TO 17 
SO TO 14 

18 IF(TREQ.LT.TOUT) SO To 17 
14 TREQ=TREQ+FNPT 
60 TO 15 
17 CONTINUE 

IF (SW(9) ) 60 TO 5*sO 
CALL CLOCK (YM) 

ZM=YM-XM 
TAB =TAB+ZM 
ISCT=N+1 

IF(NPT-N)20*21r22 

20 NCT=N-(N/NPT)*NPT 
SO TO 23 

21 NCT=0 

GO TO 23 

22 NCT=N 

23 L=K+1 

C TO COMPUTE FIRST aNU SECOND SUMS 

CALL SUMS 

IF (SW ( 12) ) GO TO 11 
WRITE(3*202) 

11 CALL CLOCK (XM) 

IF ( ISKT ( 16) ) SW(20)=.TRUE. 

52 K=L 

60 CONTINUE 
T=T+CDEL 

TOUT = T/TC+.5D-13 
ETIME=T*TC1+DSTART 
60 TO (61*62) » ITT 

61 IF(FT1+(FNP+1.D0)*CDEL .LT.T) GO TO 300 
SO TO 43 

62 IF(FT1+(FNP+1.DO)*COEL .GT.T) GO TO 300 
43 CONTINUE 

IF(ETIME.LE.DAYI) GO TO 1000 
CALL EPHQAN 
1000 CALL CSTEP 
M=K-3 

SW(12)=. FALSE. 

IF( (X(M*3) ,LT.0.0i;0.ANn.X(M-l*3) .GT.0.0D0) .AND. ( ISWT( 11) .r,R. TqWT( 
18))) CALL TNGOE 
IF (SW ( 12) ) GO TO Ol 
4l INC0WL=1NC0WL + 1 
ITERS= ITERS + ITtR 
IF ( 1TER.LE.3) 60 10 30? 
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WRITE ( 3 . 325) TOUT 

325 FORMAT (5X. 45HRUN TERMINATED. .. ITERATIONS GREATER THAN 3./ 

* 1H0 » 40X . 11HFINAL TIME= 024.12) 

GO TO 550 

3U2 IF (SW ( 19) ) IPT=IPT+1 

IFUPT .LT. bO) GO TO 303 
1PT=0 

SW(19)=. FALSE. 

M0DE=1 

IF(Ll.EQ.L2)M0DE=2 
303 IF (MODE. EQ. 4) GO TO 5f 
Sw ( b ) — • F ALSE . 

SW (9)=. FALSE. 

CALL TEST 

IF ( SW ( 9) ) GO TO 550 
IF ( SW ( 8) ) GC TO bO 
ISCT=ISCT+1 

:j 6 IF(SW(14) .AND.ISWT(ll) ) GO TO 50 
lF(ISWT(a) ) GO TO 6b 
M=MAXO(5rOR(,ER/2) 

54 TI=TREQ*TC 

GO TO ( 57 . 56) .ITT 
>;,7 IF (TI .GT . T-FNP*CDF L) Go TO 50 
oO To 59 

■)8 IF'(TI.LT.T-FNP*CJEL) GO TO 50 
59 Ml— 1 

.->3 sw ( 10 ) =. TRUE • 

N1=K 

CALL HEM1NT 
CALL OUTPUT 
TKEG— TREOi +FNPT 
GO TO 54 
o5 HCT=NCT+1 

IF(NCT-NPT)b0. 55.55 
,5 NCT=0 

CALL OUTPUT 
50 K=K+1 

IF (K .LE.200 ) GO TO 60 
00 51 1=1.3 
GO 51 K=l.NO 
U=2o0+K-N0 
X (K » I ) =X ( J » I ) 

AD(K . I )=XD( J» I ) 
si XDD (K » I ) =XDO (U. I ) 

L=NO+l 

IF (iSCT.GT.NO) ISCT=NG 
KSAVE(4 )=T-FlOAT(iSCT-1 )*CDEL 
GO TO 52 

300 Sw ( 15 ) = . TRUE , 

CALL CLOCK (YM) 

T=T-CDEL 
TM= T-OLDl 
ZM=YM-XNI 
CALL RESUME 
Sw(13)=. FALSE. 

SbO PX(N+5)=CSAVE(1) 

CX(N+3)=CSAVE(2) 

IF ( .NOT. (ISwT (11) .Ok.IS-T (18) ) ) Go TO 551 
ljO FORMAT ( lhO . 1 9HN0D/;L CROSSING DaTa/ 1H0 12X4hN0DE27X1HX27X1Hy27 x 1 HE ) 
1)1 F0RMAT(13X»I3.15x,3(5X.U24.16) ) 
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10 2 FORMAT ( 4 ( 5X r Uc4. lb) ) 

«RlTE(3r 100) 
uO 104 I=lrlNOUt 
J=1 

/»RITE(3rl01> Jr ( X.NlJUc ( I . L ) rL=lr3) 

TNOu(I)=TNOu(i)/TC 

1J4 WRlTE(3rl02) TNOO ( I ) r ( XLNODE ( I r L ) ,L=lr3) 

GO TO b51 
1 j 3 CONTINUE 
STOP 
END 

W ElT CKD1FF r 1 r 670030 r i+7b76 
ffl FjF 

SUBROUTINE CKuIFF (lu) 

c lu- integer signifying difference to be Computed. 

double PRECISION Xr xOr XDur XXDOrDlF'F r CDELr TTr T r TREO r tout 
DOUBLE PRECISION OlwCrFLrFI 

COMMON/ .VORKER/X ( 200 r 3 ) ,Xu(200r3) rXUu(200'3> rXXDD(3) HTlFF(hO»3) » 
I CUtL r TT t T r TRE 0 'TOUT iKrfJ 

DIMENSION BINC(dO) 

IF (K.GT.l) GO TO 2 
UO 1 J= 1 r 3 

1 UIFF (LD r J) = XuO ( K r J ) -XDi (K-lrJ) 

GO TO 4 

2 FL=LU 
BINC(1)=FL 
DO b J— 1 r 3 

b DIFF (LDr J) =XDO (K r j) -UIi> C ( 1 ) *XDD (k- 1 r J) 

00 3 I=2rL0 
F 1 = 1 

bINC<I)=(BINC(I-l)*(FL-FI+1.00G>)/FI 

L=K-I 

UO 3 J=lr3 

3 UIFF (LOr J) =01FF (Lu r J) + ( (-1 . 000 ) **I*BINC ( I ) *XOD(Lr J) ) 

4 RETURN 
END 

0 PUT COEFFrlr67Q830r 47576 


*0 p\jF lo) 

SUBROUTINE COLFF ( x r T r N > ORDER' A ) 


c 

c 

X 

(N) 

dependent variables 

c 

T 

(N> 

INDEPENDENT VARIABLES 

c 

N 


NUMBER of INPUT VARIaBLfS (MAXIMUM OF 2o> 

c 

Or(UER 


ORDER uF POLYNOMIAL FIT (MAXIMUM OF 7) 

c 

c 

A* 

(ORDER+1) 

COEFFICIENTS OF POLYNOMIAL FIT 


bOUbLE 

PRECISION 

X 1 1 ) 


OOUbLE 

PRECISION 

r(l> rA(l) rSUMrBTrBTB 


COMMON/SCRTCH/BT(«rEO) ,BTB(8'8) rSUMrMrGl (49) 
INTEGER ORDER 
iM=MIN0 ( ORDER 41 r 6) 

UO 100 1=1 fN 

BT(lrI)=1.0D0 
DO 100 0=2 r M 

100 BT(JrI)=T(I) *BT ( J-l r 1 ) 

UO 300 1=1 r M 
DO 300 J=1 r M 
SUM=0 • 0D0 
DO 200 K=lrN 

200 SUM=SUM+BT(I,K)*BI (J'K) 
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300 bTb( I , J)=SUf" 

CALL DNVERT(M»BTJ,8,A) 

DO 500 1=1, M 
A(1)=0.0D0 
DO bOO J=1,N 
SUM=0.000 
DO 400 K=1,M 

40C SUM=SUM+BTB(1,K)*hTU\,b> 

50C A(l)=A(I)+SUi-1*X( J) 

RETURN 

END 

» Ei_T CSTEP, 1,670830, 47bb0 
0 EOF i«i 

SUBROUTINE CSTEP 

DOUBLE PRECISION X,XO,xDO»XXbD,DIFF»COEL»TT»T»TREQ»TOUT 
DOUBLE PRECISION Sl,S2,PX,PXD»CX,CX0,CT0L,SAVE 
DOUBLE PRECISION &Um 1»SUM2 
DOUBLE PRECISION sC0,SPCb 

COMMON/WORKEK/X (200,3) , Xu (200, 3) ,XDO( 200,3) ,XX0D(3) » NIFF (60,3) , 
l CUEL,TT,t,TR£u,TOUT»K,N 

COMMON/COWS/SI (3) ,S2(3) ,PX(63) ,PXD(63) ,CX (63) *CXD(63> ,CT01 , SAVE (30 
l ,3), ITER 

COMMON/TIMES/ TAB, COWL, FORCS 
DIMENSION SC 0 ( 3 ) , SPC 0(3) 

INTEGER B 
CALL CLOCK (Xt_) 

DO 1 1=1,3 
SUM1=0.D0 
SUM2=0.D0 
DO 2 U=1,N 
B=N-J+1 

SUM1=SUM1+PX ( d+3 ) *D IFF ( 8 * I ) 

2 SUM2=SUM2+PXU (B+2) ♦DlFF (B , I ) 

SCO ( 1 ) =SUM1 + PX ( 3 ) * XhD(K-lrl) 

X (K, I ) = SCO(I) + S2(I) 

1 XD(K,I)= (SUM2 + PXD(2)*XDD(R-l,I)+Sll I) )/CDEL 
CALL FRCS 
DO 15 1=1,3 
SAVE ( 1 » I ) =DIFF ( 1 * I ) 
ulFF(l,I)=XDU(K,I) - Xi,D(K-l,I) 

DO 7 J=2,N 
SAVEtU, I ) =OIFF <U, I ) 

7 DIFF(J,I)=DIFF(J-1,1)-SaVE(J-1,I) 

15 CONTINUE 
1TER=0 

100 DO b 1=1,3 

SPC0(I)=SC0(I) 

SUM1=0.D0 
SUM2=O.DO 
DO 3 J=1,N 
B=N-J+1 

SUM1=SUM1+CX (B+3)*DIFF (B, I ) 

3 SUM2=SUM2+CX0(B+2)*UlFF(B»I) 

SC0(I)=SUM1+CX(3)*XDD(k,I) 

X(K,I)=SC0(I) + S2(l) 

5 XD(K,I)= <SUM2-CXD(2>*xDU(K,l) + Sl(I))/COEL 
1TER=ITER+1 
DO 10 1=1,3 

IF(uABS(SPCOm-ScO(I) ) .LE.CTOL) GO TO 11 
10 CONTINUE 
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CALL FRCS 
DO 13 1=1*3 

DIFFU,I)=XDQ(K*I)-XDD(K-1 .1) 
DO 13 J=2 * N 


13 DIFF<J*I)=DIFF(J- 

l.D- SAVE(U— 1*1) 



IF (ITER .LE. 3) 
RETURN 

GO TO 100 



11 DO 12 

1=1.3 




S1(I) 

=S1(I) + XDO(K, I ) 



12 S2U> 

=S1 ( I ) +S2( I ) 




CALL 

CLOCK (XR) 




COWL= 

COWL + XR - 

XL 



RETURN 




END 





CUT OIFF » 

1,670830* 47h81 



EOF H 




e 




DIFF 

7 

PURPOSE 

TO CALCULATE THE DIFFERENCE BETWEEN ANY TWO TIME 

DIFF 

8 


POINTS IN THE 20TH CENTURY 

DIFF 

9 




DIFF 

10 




DIFF 

11 

CALLING SEQUENCE 

CALL DIFF<IYMD1*IHMS1*IYMD2»IHMS2»I0AY»tSEC> 






OIFF 

13 

SYMBOL 

TYPt 

DESCRIPTION 

DIFF 

14 




DIFF 

15 

IYMDK1) 

I 

input - date in the form yymmdo 

OIFF 

17 

ihmskd 

I 

INPUT - TIME ON DATE IYMDl IN FORM YYmMOD 

DIFF 

19 

IYM02U) 

I 

input - date in the Form yymmdo 

DIFF 

21 

IhMs2(1) 

I 

INPUT - TIME ON DATE 1YM02 IN FORM HHMMSS 

DIFF 

23 

luAY(i) 

I 

OUTPUT - ELAPSED FULL DAYS DIFFERENCE 

DIFF 

24 



IDAY IS NEGATIVE IF IYMD2»IHmS2 
IS THE EARLIER TIME 

DIFF 

26 

ISEC(l) 

I 

OUTPUT - REMAINDER OF DIFFERENCE IN SFCONnS. 

DIFF 

29 



ISEC FOLLOWS THE SIGN OF IDAY 

DIFF 

32 




DIFF 

33 

ROUTINES 

REQUIRED 

RYiVOI 

DIFF 

39 


SUBROU TINE 0 IFF ( I Y Mu 1 » 1 HMS 1 » I YMD2 » I HMS2 * I DAY » ISEC ) 


DIMENSION NYEAR(12) *LYt- AR(12) 

SET THE ELAPSED DAYS AT THE BEGINNING OF EACH MONTH 

SINCE THE START OF A NON-LEAP YEAR 

DATA NYEAR /0. 31*59*90, 120. 151*181. 212*243*273*304, 334/ 

set the elapsed days at the beginning of each month 

SINCE THE START OF A LEAP YEAR 

DATA LYEAR /0* 31*b0 .91 , 121 . 152, 182*213,244,274*305,335/ 

IYEAR1=0 

1YEAK2=0 
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o O O OOO OOO OOO OOO OOO OOO OOO O D 


check for « diffei Ei.ce of less than one day 

IF 1 1 YMD1.EQ. 1 Y MU2 ) GOT( 4000 

SEPARAT 1 IYl-ibl AN, IYMD2 INTO THREE i.OkuS EACH 

CALL RYMOI ( 1 YMD1 » 1 Yi # I...1 * ID1 ) 

CALL RYMUl (IYMD2* IY2» I:v ?» I 02) 

COMPUTE THE ELaPSM, DAYS SINCE JANUARY 0»1900 

IYEaR 1=3652B*(IY1-1)/1; j U 
IYEAK2=3b525*lIY2-l)/l ,0 

CORRECT THE VALUES FOR LEAP YEAR 

IF (MOU(IYlfZ) .EU.0) oOTO 1000 
IYEAR1 =IYEAIU+NYEaR( I.-lD + IDl 
5C J IF (NiODt IY2»4) .ECi.o) oOTO 2000 
IYEaR2=IYEAR2+NYL/vR(1m. >+ID2 
OOTu 3U0u 

1000 1 YEARl = I YEAP x+LYEaR l Ini ) +101 
UOTo bOO 

20 jO I YEmR2=I YEAR2+LYEAR ( IMP ) + ID2 

convert llapslu d.:Ys into elapsed seconds 

3000 IYE«Kl=l YEARl*«o400 
lYEAR2=IYEAF2*6b400 

CALCULATE ELAPSED SECONDS INTO EACH DAY 

4000 ISECl=IHMSl-40*(lhMsl/l00)-2400*( IHMSl/10o00> 
ISEC2=IHMS2-4g* ( I mMs2/10C ) -2400* ( IHMS2/10000 ) 

subtract the uo t lapsed seconds values 

ISEC=I YEAR2+1SEC2-I YlA>. 1-ISEC 1 

COMPUTE IOaY 

IDAY=ISEC/86200 

COMPUTE ISlC 

ISEC=ISEC-lDrt Y*tib400 
C 

RETURN 

C 

END 

Q FlT DNVERT » 1 » 670830 » 47dB2 
8 EOF >'«» 

SUBROUTINE ONVER r ( Nr A» ijT » I ROW) 

UOUbLE PRECISION A>X»Y,Z 
DIMENSION A(NT»NT) »lROi. (NT) 

DO 1 I = 1 * N 

1 IROw(I) = I 

L = 1 

2 IF(L.EO.N) DO TO b 
I = L 
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K = L + 1 
DO 4 J - KrM 
X = ABS(A(J,L>) 

Y = ABS(A(I»l) ) 

IF ( Y .GT.X) GO TO 4 

1 = J 

4 CONTINUE 
IFU.EQ.L) GO TO b 
DO 8 J = 1 f N 

2 = A(L* J) 

A ( L » J) = A ( I f J) 

8 A ( 1 » U) — 2 
K = I ROW (L) 

IROw(L) = iKOto(I) 

IROw(I) = K 

5 IF (A(L'L>-0.) 9»b»9 
b PRINT 10 

10 FORMAT ( 19H MATRIX IS SINGULAR) 
RETURN 

9 2= 1.0 / A(LM-) 

IF(L.LE.l) GO TO 13 
K = L-l 

DO 11 I = 1 * K 

11 A(L.I) = A(i_ri) * 2 
IF (N.LE.L) GO Td )5 

13 K. = L + 1 

DO 14 I = K »N 

14 A ( L * I ) - A(Lrl) * 2 

15 IF(L.GE.U) GOTO 23 

DO 22 I = K » N 
IF (A(1»l)-0.) 16»22» 1h 

18 IF(L.LE.l) GOTO 20 
DO = L-l 

DO 19 J = 1 t o J 

19 A ( 1 r J) = A ( 1 r J) — a ( 1 r L ) * A(L»j) 

20 JO = L + 1 

DO 21 J = JJ* N 

21 A(I»J) = A ( l » J) - A(I»L) * A(L*J) 
A ( I r L) = - A ( I ’ L) * 2 

22 CONTINUE 

23 A(L»L) = 2 

IF (n.LE.L) ou TO 2b 
L — L + 1 
GO TO 2 

ib DO 32 L -■ 2 > U 
K = L - 1 
DO 02 I = 1»K 
IF<A<I»L)-0.) 2b»o2»2b 
2b DO 27 J=1*K 

27 A(I.J) = A(Ifd) - A ( I » l ) * A ( L ► J ) 
IF(N.EO.L) GOTO 31 
i\K = L + 1 

DO 00 J = KK.»N 

AO A ( I f J) = A ( I * J) - A ( I » 1. ) * A (L » J) 

.11 A ( I » L) = - A II*L> * A(L»L) 

12 CONTINUE 
DD — N — 1 

DO 37 I = 1 rJJ 

IF ( I ROW ( I ) .i 0 • I ) GO TO 37 

13 DO 34 J = I»u 
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3. c 


IF <IR0*(U) .FU.l) GG 10 

)4 CONTINUE 

,lb [ill J6 K ^ liN 

Z = A (K* I ) 

A ( K * I ) = A ( K * U ) 

3b A(K*J) = Z 

IKCw ( J) — IKOa( I) 
t>l CONTINUE 
HETURN 
END 

a f I_T DRAG* 1 * 871807 1 48-, 6u 

a e of ioj 

SUDkOUTINE DKrtC 

COMMON/ WORKEK/X (2b0 * 3) ,*u<200*3> *xn0(2)0*3) *XXDD<3) *DIFF (nfl * 3) * 
* CDEl * |T*T*IRE«*ToUT*K*N 


COMMON/CONST id/GM# At * R * l?So * RO * THETg . TCI 
C0MM0N/C0NST3/EMASS(2) , XYZ14) *LS 

COMMON/PERTS/ALT ( 4 1 ) *0fN(41*2) *H..,(40*2) » rtE * ESQ * >J * CN* COP * RmOI . 


OOUuLE 

DOUbLE 

OOUbLE 

4c 

* 


rtH02 * Kh03 * Rh04 * C APR * ARE A * SATmAS * CD * MU 
PRECISION x *XD* xDU* XXOO* D IFF* COEL* TT* T* TrEO* TOUT 
precision gm* ae»k* RSG *Ru*THETG*TC1*EMASS*XYZ 

PRECISION ALT * DEN* DRAG ( 3 ) * VR * wE* 03 * Q4 * G5 * ESQulEL * R7 * B * FTA j 
liI*HM,Rh01*KH02.RH03*RH04*C0SPSI*RH0*RMA**RMlN* 
CN ► COP *CAPR* AREA* SATMAS* CO *DEXP 


DOUbLE PRECISION OSuRT.UC0S*DSlN 
DRAG(1)=XD(K* 1)+WE*X(K*2) 

DRAG (2) =XD (K * 2 ) -WE*X (K * 1 ) 
URA6(3)=XD(K*3) 

VR=0.00 
DO 4 1 = 1*3 
4 VK=VR+0RAG(I)**2 


VR=DSGRT(VR) 

G3=X (K* 1 ) **2+X ( K * 2) **2 

G4=Q3/RSG 

G5=ESG*Q4 

DEL=(ESQ*X(K*3)*USQRT(Q3> ) / (RSG-ESQ*Q3) 
K7= B*(l.DO+O.5D0*Q5*(l.D0+O.75DO*Q5) ) 
ETA=QEL*R7/R 

HI=(R-R7)*(1.00-ETA*DEL/2.0DU)*CAPR 
IF(HI.LT.1.0D+03) GO TO 22 
Rh0=0 .DO 


GO TO 21 
22 DO 10 1=1*41 

IF (Hl.GT.ALT(I) ) GO TO 10 

RMAX= DEN (1-1*1) *DEXP ( (ALT (I-1)-HI)/HM (1-1*1) ) 

RMIN= DEM 1-1*2) *DEXP( (ALT ( I -1) -HI )/HM( 1-1*2) ) 

60 TO 15 
10 CONTINUE 

15 C0SPSI=X(K*1)/R*(XYZ(1)*DC0S<RH04) ) 

* +X(K*2)/R*(XYZ(2)*DSIN( RH04 ) ) 

* +X(K*3)/R*XYZ(3) 

COSPSI=(USQRT( ( 1 .DO+COSPSI ) /2 . DO ) )**MD 
RHO=(RMlN+(RMAX-RMlN) *COSPSI ) *1 . 0D-10*CAPR 

21 DO 20 1=1*3 

20 DRAG ( I ) = ( 1 . DO+RHOl ) *CDP* ( 1 . DO+RH02+T ) *RHO*VR*DRAG ( I ) 


DO 25 1=1*3 

25 XXDD ( I ) =XXDD ( I ) -DRAG ( I ) 
RETURN 
END 

fi ELT EGRAV *1*670830* 47574 
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a EOF u 

SUBROUTINE EGRAV 
C 

c****notatiqn**** 

c 

C PSI-LATITUDE (GEOCENTRIC) 

C LAMBDA-LONGITUDE (+EASTWARD) 

C R-GEOCENTRIC RADIUS TO SATELLITE IN EARTH RADII 
C GM-GRAVITATIONAL CONSTANT TIMES MASS OF EARTH 
C P(M.N) -COEFFICIENTS OF LEGENDRE POLYNOMIAL 

C C (N.M) “COEFFICIENTS OF COSINE FUNCTION 
C S(N»M)-COFFICIENTS OF SINE FUNCTION 

C INDEX 1-DEGREE OF SUMMATION PLUS 1 

C 

DIMENSION P(22.20) .S(20.23) . COSLAM ( 21 ) , SINLAM < 21 ) .TPSIM(2l) » 

. C (20.23) 

C 

DOUBLE PRECISION X ► XD » XDD. XXDD.DIFF.CDEL . TT. T» TREQ.TOUT 
DOUBLE PRECISION GM » AE ► R » RSQ » RQ . THETG * CS . S » C . R ASAT » PR . DR * PLAMDA » 

* DLAMDA .PPSI fOPSI » LAMBDA . GMR . PRR . PLX Y . PPTP . COSLA 

* SINLAM. TPSIM.TANPSI .SINPSI .COSPSI » ZERO. CP3.CL2. 

* F1.F2.F3.F4.FN1.RN.RINV.FM .PI. XYSQ.RTXYSQ»DX(3) 

* »DATAN2.DS1N»DC0S.DSQRT *P»TC1 
C 

COMMON/ WORKER/X (200.3) » XQ ( 200. 3) .XDD(200.3) .XXDD(3) »DIFF(60.3) . 

* cuel.tt.t.treq.tout.k.noro 

C0MM0N/C0NST2/GM.AE.R .RSQ. RQ» THETG. TCI 
C C AND S COEFFICIENTS ARE STORED IN SAME MATRIX 
COMMON/FMODEL/CS (20.23) . INDEX1 , INDEX3 
C 

EQUIVALENCE (TPSIM(2) .TANPSI) . (C.S.CS) » (P.SINPSI) . (P(2»l) .COSPSI), 
. (TPSIM.ZtRO) 

DATA ZERO * SINLAM ( 1 ) .COSLAM(1)/2*O.DO»1.DO/ 

C SET ZERO ELEMENTS OF LEGENDRE POLYNOMIALS 

DATA P (3* 1) *P (4. 2) * P( 5. 3) * P (6*4 ) »P(7» 5) * P ( 8.6) .P (9.7) * 

. P(10.8).P(11.9)*P(12.10).P(13.11)*P(14.12)*P(15»13) .P(16.14), 

. P(17.15) .P( 18.16) ,P( 19.17) .P (20.18) .P (21.19) .P( 22. 20)/ 

. 20*0. DO/ 

C 

RASAT=DATAN2(X(K.2) .X(K.i) > 

XYSQ=X(K.1)**2+X(K.2)**2 
RTXYSQ=DSQRT(XYSQ) 

LAMBDA=RASAT - THETG 
SINLAM ( 2) = DSIN( LAMBDA) 

COSLAM ( 2 ) =DCOS ( LAMBDA ) 

SINE. COSINE .AND TANGENT OF LATITUDE 

SINPSI=X(K.3)/R 
COSPSI=RTXYSQ/R 
TANPSI= SINPSI/COSPSI 

RINV=AE/R 


****ALL INDEXES 1 HIGHER THAN IN ANALYSIS**** 


calculate polynomial terms 

...NOTE - P TAKES form P(M.N) 


INDEX2=INDEX1-1 
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CP3=3.UO*COSPSI 

P(1*2)=1.5D0*SINPSI**2-.5Q0 

P(2*2)=CP3*SINPbI 

P(3*2)=CP3*C0SPSI 

TPSIM<3)=2.D0*TANPSI 

CALCULATE AND SAVE SINES AND COSINES OF LONGITUDE 

CL2=2.O0*COSLAM(2) 

SINLAM(3)=CL2*SINLAM(2) 

COSLAM ( 3) =CL2*C0SLAM ( 2 ) -1 . DO 

DO 120 N=3* INUEX2 

F1=N 

F2=F1-1.00 

F3=2.DO*F1-1.UO 

F4=F3*C0SPSI 

N1=N-1 

N2=N-2 


ZONAL HARMONICS (M = 0) 

P ( 1 * N ) = (F3*SINPSI*P (l*Nl)“F2*P(l*N2) )/Fl 
IF ( 1N0EX3.LT >2) GO TO 120 
NX=MINO(N*INO£X3) 

DO 110 M=2»NX 

TESSERAL HARMONICS (M NON ZERO* LESS THAN N) 

110 P(M*N)=P(M*N2)+F4*P(M-1*N1) 

IF(NX.LT.N)GO TO 120 
NN1=N+1 

SECTORAL HARMONICS (M EQUAL TO N* NON ZERO) 

PINN1*N)=F4*P(N*N1) 

TPSIM(NN1)=TPSIM(N)+TANPSI 

SINLAM ( NN1 ) = CL2*SINLAM ( N) -SINLAM (N1 ) 

C0SLAM(NN1)=CL2*C0SLAM(N) -COSLAM (Nl) 

120 CONTINUE 

INITIALIZATION FOR SUMMATION FOR PARTI ALS 


140 PR=0. 

PLAMDA=0. 

PPSI=0. 

FN1=2. 

RN=RINV 

SUMMATION FOR PARTI ALS 

DO 2S0 NC=2 * INDEX2 

NS=21-NC 

RN=RN*RINV 

FN1=FN1+1.D0 

FM=0. 

DR=0. 

DLAMDA=0. 

DPSI=0. 

N1=MIN0(NC+1*INDEX3) 
DO 200 MC=1*N1 
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MS=24-MC 

Pl=P(MCrNC) 

IF(MC.EQ.1)G0 TO 150 
FM=FM+1.00 

PARTIAL WITH RESPECT TO LAMBDA (SUMMATIOM) 

DLAMDA=DLAMDA+FM*Pl* ( S ( NS t MS ) *C0SLAM ( MC ) -C ( NC t MC ) *SINLAM ( MC ) ) 
150 F 1=C (NCrMC) *C0SLAM ( MC ) +S ( NS t MS ) *SINLAM ( MC ) 

IF(Fl)175r200»175 

PARTIAL WITH RESPECT TO R (SUMMATION) 

175 DR=DR+F1*P1 

PARTIAL WITH RESPECT TO PSI (SUMMATION) 

DPSI=DPSI+Fl*(P(MC+lnNC)-TPSIM<MC)*Pl) 

200 CONTINUE 

PK=PR+DR*FN1*RN 
PLAMDA=PLAMDA+DLAMDA*RN 
250 PPSI=PPSI+DPSI*RN 
GMR = GM/R 

COMPLETE PARTIAL WITH RESPECT TO R 
PR=-GMR* ( 1 . 0D0+PR ) /R 
COMPLETE PARTIAL WITH RESPECT TO LAMBDA 
PLAMDA=GMR*PLAMDA 

COMPLETE PARTIAL WITH RESPECT TO PSI 
PPSI=GMR*PPSI 

CONVERT ACCELERATION IN SPHERICAL COORDINATES TO ACCELERATION IN 
RECTANGULAR COORDINATES (MULTIPLY BY MATRIX OF PARTI ALS OF 
SPHERICAL WITH RESPECT TO RECTANGULAR) 

PRR=PR/R 

PLXY=PLAMDA/XYSO 

PPTP=PRR-PPSI*X(K»3)/(RTXYSQ*RSQ) 

DX ( 1 ) =X ( K * 1 ) *PPTP-PLXY*X ( K 1 2 ) 

DX(2)=X(K»2)*PPTP+PLXY*X(Krl) 

DX(3)=PRR*X(K»3)+PPSI*RTXYSQ/RSQ 
XXDD ( 1 ) =DX ( 1 ) 

XXDD(2)=DX(2) 

XXDD(3)=DX(3) 

RETURN 

END' 

0 ELT EPHtMr 1 ► 670927 » 50403 
IS EOF U 


c 

c 

CALLING 

SEQUENCE 

CALL EPHEM ( I Yl » D1 * AO) 



c 

I Yi 

I 

INPUT 

YEARS SINCE 1960 



c 

D1 

□ 

INPUT 

TIME OF INTEREST IN DAYS 

SINCE-IY1- 

c 

AO 

R 

OUTPUT 

VECTOR OF NINE ELEMENTS 



c 




ELEMENTS 1-3 INERTIAL 

Xf YrZ 

COMPONENTS OF 

c 




UNIT VECTOR TO 

MOON (METERS) 
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ELEMENT 4 

ELEMENTS 5-7 

ELEMENTS 8 
ELEMENT 9 


RANGE TO MMON (METERS) 
INERTIAL X#Y#Z * COMPONENTS OF 
UNIT VECTOR TO SUN (METERS) 
RANGE TO SUN (METERS) 
EQUATION OF THE EQUINOXES 


SUBROUTINE EPHEM ( I Y1 * 01 * AO ) 

DOUBLE PRECISION D * T * TWOPID * CRAD*DRSEC * THDOT 1 * THD0T2* DUMMY 
DOUBLE PRECISION D1*DMIN*ETIME 
DOUBLE PRECISION SLPT 
OIMENSION AO (9) 

DIMENSION DAYS(IO) 

COMMON/CONST /DRAD* TWOPID * DRSEC * RAD * RSEC 
C0MM0N/C0NST1/DUMMY ( 10 ) * THDOT I * THD0T2 * OMIN * ETIME * I YBEG 
COMMON/CSUN/ASUN # PMOON # EMOON * T ASUN (10) 

REAL LS*LM 

EQUIVALENCE (SL*SLMGM) * (SP.SLSGS) * (CP*CLSGS) 


DATA DAYS/O. *366. *731. #1096. *1461. #1827. *2192. #2557. *2922. *3288. 
tPHEMERIS DAYS SINCE JAN 0.5 1900 
D=D1+21913. 5D0+DAYS ( I Y1 ) 

TIME IN JULIAN CENTURIES SINCE JAN 0.5 1900 

T=D/36525.0D0 

TSQ=T**2 

DSG=1.0E-8*D**2 

MEAN OBLIQUITY OF THE ECLIPTIC 

EPS=RAD* (23. 452294-. 0130125*T-. 164E-5*TSQ ) 

MEAN LONGITUDE OF THE SUN 

LS=DMOD (DRSEC* ( • 10 0690 80407+ . 1296027681 3D9*T+1 • 089*TSQ) * TWOPID) 

SLS=SIN(LS) 

CLS=COS(LS) 

S2LS=2.*SLS*CLS 

C1=CLS**2 

C2LS=2.*C1-1. 

S3LS=SLS*(3.-4.*SLS**2) 

C3LS=CLS*(4.*Cl-3. ) 

MEAN LONGITUDE OF PERIGEE OF THE SUN 

GS=RSEC*( 1012395. 0D0+6189. 03* T+1.63*TSQ) 

SGS=SIN(GS) 

CGS=COS(GS) 

MEAN LONGITUDE OF THE MOON 

LM=DMOD( DRAD* (270. 434164+13. 1763965268D0*D-.85E-4*DSQ) .TWOPID) 


54 



ooo ooo ooo o ooo o ooo 


SLM=SIN(LM) 

CLM=COS(LM) 

S2LM=2 • *SLM*CLM 

C1=CLM**2 

C2LM=2.*C1-1. 

S3LM=SLM* (3.-4. *SLM**2 ) 

C3LM=CLM*(4.*Cl-3. ) 

MEAN LONGITUDE OF LUNAR PERIGEE 

GM=DMOD (DRAD* ( 334 . 329556+ . 1 114040803D0+D-. 7739E-3*DSQ ) . TWOPID ) 

SGM=SIN(GM) 

CGM=COS(GM) 

LONGITUDE OF MEAN ASCENDING NODE OF THE LUNAR ORBIT 

OM=DMOD (DR AD* (259. 183275-. 529539222E-1*D+ • 1557E-3*0SQ) .TWOPID) 

S0M=SIN(0M) 

COM=COS(OM) 

S20M=2.*S0M*C0M 

C20M=2.*C0M**2-1. 

C1=SGS*CLS 

SLSGS=SLS*CGS 

SLSPGS=SLSGS+C1 

SLSGS=SLSGS-C1 

CLSGS=CLS*CGS+SLS*SGS 

SLMGM=SLM*CGM 

C1=SGM*CLM 

SLMPGM=SLMGM+C1 

SLMGM=SLMGM-C1 

CLMPGM=CLM*CGM+SLM*SGM 

S2LM0M=S2LM*C0M-S0M*C2LM 

S2LS0M=S2LS*C0M-S0M*C2LS 

C2LM0M=C2LM*C0M+S2LM*S0M 

S3LMGM=S3LM*CGM-SGM*C3LM 

S3LMGS=S3LM*CGS-SGS*C3LM 

S3LSGS=S3LS*CGS-SGS*C3LS 

C3LSGS=C3LS*CGS-S3LS*SGS 

S2LSGL=S2LS*CLMPGM-SLMPGM*C2LS 

NUTATION IN LONGITUDE 

DELPSI=DRSEC*( ,2088*S20M-SOM* ( 17.2327+. 01737*T) -1 .273*S2LS 
-• 2037*S2LM+. 1259*SLSGS-. 0496*S3LSGS+. 0214+SLSPGS 
+. 0675+SLMGM-. 0342*S2LMOM-. 0261+S3LMGM+ . 0114*SLMPGM 
+.0124*S2LS0M-.0149*S2LSGL) 

NUTATION IN OBLIQUITY 

DELEPS=DRSEC* ( 9.2106+COM-. 0904*C20M+ • 552+C2LS+. 0884*C2LM 
+.0183*C2LMOM+.0216*C3LSGS) 

TRUE OBLIQUITY OF THE ECLIPTIC 

AEPS=EPS+DELEPS 

SOB=SIN(AEPS) 

C0B=SQRT(1.-S06**2) 
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EQUATION OF EQUINOXES 

AO(9)=COB*DELPSI 

ECCENTRICITY OF EARTH’S ORBIT 

ESUN=.01675104-.418E-4*T 

APPARENT LONGITUDE OF THE SUN 

03=DINT(Ul) 

D2=U1-D3 

D2=D3*TH00T1+l>2*THD0T2 

ALS=LS+2.*ESUN*SIN(U2-TASUN< IY1) ) 

SALS=SIN(ALS) 

CALS=COS(ALS) 

INERTIAL UNIT VECTOR TO SUN 

A0(5)=CALS 
AO(6)=SALS*COB 
AO(7)=SALS*SOB 

RAOIUS VECTOR TO SUN 

PSUN=ASUN* ( 1 . -ESUN**2 ) 

SLPT =CALS*CGS+SALS*SGS 
AO ( 8 ) =PSUN/ ( l . +£SUN*SLPT ) 

TIME IN EPHEMERIS DAYS SINCE JAN 0.5 1960 

T1=D-21914. 

C 

CL=CLM*CGM+SLM*SGM 

SF=SLM*COM-SOM*CLM 

CF=CLM*COM+SLM*SOM 

SD=SLM*CLS-SLS*CLM 

CO=CLS*CLM+SLS*SLM 

S2L=2.*SL*CL 

C2L=2.*CL**2-1. 

S2D=2. *SU*CD 
C2D=2.*C0**2-1. 

S2F=2.*SF*CF 

C2F=2.*CF**2-1. 

SLPP=SL*CP 

C1=SP*CL 

SLMP=SLPP-C1 

SLPP=SLPP+C1 

SLP2D=SL*C2D 

C1=S2D*CL 

SLM20=SLP2D-Ci 

SLP20=SLP2D+C1 

CLP20=CL*C2D 

C1=SL*S2D 

CLM2D=CLP2D+C1 

CLP2D=CLP2D-C1 

SPP2D=SP*C2D 

C1=S2D*CP 
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SPM2D=SPP2D-C1 
SPP2U=SPP2D+C1 
C1=SL*C2F 
C2=S2F*CL 
S4D=2.*S2D*C2D 
C4D=2.*C2D**2-1. 

S=LM-OM+.1096*SL-.0222*SLM2D+.0115*S2D+.0037*S2L 

APPARENT LONGITUDE OF MOON 

ALM = 206265. *LM+22640.*SL-4b86.*SLM2D 
. +2370.*S2D+769.*S2L-668.*SP-412.*S2F 

. -212.*(S2L*C2D-S2D*C2L)-206.*(SLPP*C2D-S2D*(CL*CP-SL*SP) ) 

. +192.*SLP2D-165.*SPM2D+148.*SLMP-125.*SD 

ALM =(ALM-109.*SLPP-55.*(S2F*C2D-S2D*C2F)-45.*(C1+C2) 

. +40.*(C1-C2)-38.*(SL*C4D-S4D*CL)+144.*( . 75*SL-SL**3> 

. -62.*SLM2D*CLM2D+28.*(SLM2D*CP-SP*CLM2D> 

. -24.*SPP2D+19.*(SL*CD-SD*CL)+18.*(SP*CD+SD*CP) )/206265. 

SALM=SIN(ALM) 

CALM=COS(ALM) 

APPARENT LATITUDE OF MOON 

ALAT-(1852 0.*SIN(S)*(1.-.00293*SIN (1.40380 8-, 000 9242203+T1) ) 

. -31 . * ( SF*CLP2D-SLP2D*CF ) -25 . * ( SF*C2L-S2L*CF ) 

. -23 . * ( SPM2D*CF+SF* ( CP*C2D+SP*S2D ) ) +21 . * ( SF*CL-SL*CF ) 

. -526.* (SF*C2D-S2D*CF) +44.* (SLM2D*CF+SF*CLM2D) )/206265. 

SALAT=SIN(ALAT) 

CALAT=SQRT ( 1 ,-SALAT**2) 

INERTIAL UNIT VECTOR TO MOON 

Cl=CALAT*SALM 
A0(1)=CALAT*CALM 
AO(2)=C1*COB-SALAT*SOB 
AO(3)=C1*SOB+SALAT*COB 

RADIUS VECTOR TO MOON 

A0 ( 4 ) =PMOON/ ( 1 . +EMOON*CL ) 

RETURN 
END 

IS ELT EXTCF. 1 »t>70830 » 47586 
« EOF IS 

SUBROUTINE EXTCF 

DOUBLE PRECISION XNODE . XDNODE * TNOD 

DOUBLE PRECISION B»C . ALPHA .FN> FI . SUM. SUM1 . CC 

DIMENSION B(20) 

COMMON /NODE/ XNODE (200.3) * XDNODE (200.3) .TNOD(200) .C(20) rCC(20) . 
IRK » J1 » ND1FF. NEXT »KEXT r INODE 
N=NEXT 
K=KEXr 
FN=N 

ALPHA=-1.D0/FN 
B ( 1 ) - ( 1 .DO- ALPHA ) /2. DO 
C(l)= l.DO-B(l) 

CC(1)=-B(1) 

DO 1 1-2. K 
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FI=I 

B(I)=(FI-ALPHA)*B(I-l)/(FI+l.O0) 

SUM=O.D0 

SUM1=0.00 

L=I-1 

BO 2 J=1»L 
M=I-J 

SUM1=SUM1+B ( J ) *CC ( M ) 

2 SUM=SUM*B<J)*C(M> 
CC(I)=-(SUMX+B(I) ) 

1 C(I)=l.DO-(SUM+b(I>> 

RETURN 

END 

S ELT FIT. 1*670830 . 47587 
0 EOF «i 

SUBROUTINE FIT(N*0R0ER1 *A*T .X) 


N 

OROLK1 

A 

(0RDER1.N) 

NUMBER OF POLYNOMIALS TO FIT 
ORDER PLUS ONE OF POLYNOMIALS 
COEFFICIENTS OF POLYNOMIALS 

T 

X* 

(N) 

POINT AT WHICH TO EVALUATE POLYNOMIALS 
VALUE OF POLYNOMIALS AT *T* 

INTEGER ORDER 1 



DOUBLE PRECISION A(5.9> .T1*X(9) .T 

Tin. ooo 

00 100 J=1»N 
100 X(U)=A(1.J> 

00 200 1=2. ORDER 1 

T1=T*T1 

00 200 J=1.N 

200 X( J)=X(J)+T1*A(I. J) 

RETURN 

ENO 

9 ELT FKOIFF. 1.670830. 47588 

a eof a 

SUBROUTINE FKD IFF ( ARRAY .OIFF.K.M.N) 

DOUBLE PRECISION ARRAY. DIFF.FK.BINC. FI .SUM 
DIMENSION ARRAY (100. 3) »DIFF(20.3) .BINC(20) 

IF(K-2>1»2.2 

1 DIFF(K.N)=ARRAY (M*M -ARRAY (M-l.NJ 
60 TO 4 

2 FK=R 

BINC(1)=FK 

SUM=ARR A Y < M . N > -B I NC 1 1 > * ARR A Y t M-l . N ) 

DO 3 1=2. K 
FI=I 

BINC(I)=(BINC(I-1)*(FK-FI+1.000))/FI 

L=M-I 

3 SUM=SUM+ ( (-1. UDO) **I*BINC ( I ) *ARRAY (L.N) ) 

OIFF (K*N)=SUM 

4 RETURN 
ENO 

8 ELT FRCS. 1.671009. 48541 

a eof a 

SUBROUTINE FRCS 

DOUBLE PRECISION X . XD . XDD * XXUD . D IFF . CDEL . TT * T * TREQ . T OUT . THETGO . 

* TH00T1 . THD0T2 • DMIN . ETIME . 6M . AE . R . RSQ . RQ . THETG . 

* C . OA Y1 . CENTER . DA Y2 . VARD ( 9 ) . EQ . TTRANS . DA Y . FDAY . 

* DSTART.DSQRT.TC1.T1 
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DOUBLE PRECISION T0L1.T0L2.TC 
DOUBLE PRECISION XYZ.EMASS 

COMMON/WORKER/X (200.3) .XD(200.3) .XDD(200.3) .XXDD(3) »DIFF(60.3) » 

* CDEL*TT .T .TREQ.TOUT .K.NORD 

C0MM0N/LIMITS/T0L1.T0L2.TC. ISCT. ISWT(30) »SW(20) . ORDER » LI »L2. MODE 
COMMON/CONST 1/THETGO ( 10 ) . THDOTl * THD0T2. DMIN.ETIME. I YBEG 
C0MM0N/C0NST2/GM . AE . R . RSQ » RQ » THETG ► TCI 
COMMON/ C0NST3/EMASS ( 2 ) ,XYZ(4) *LS 
COMMON/COFIT/C (5, 9) »DAY1. CENTER .DSTART 
COMMON/TIMES/TAB. COWL. FORCS 
LOGICAL ISWT . SW 
EQUIVALENCE ( VARD (9) »EQ) 

CALL CLOCK (XL) 

RSQ=0.D0 
DO 5 1=1.3 
5 RSQ=RSQ+X(K.I)**2 
R=DSQRT(RSQ) 

RQ=RSQ*R 

IF ( ISWT (21) ) GO TO 10 
IYl=IYBEG-59 
DAY2=DSTART+T*TC1 
TTRANS=DAY2-CENTER 
T1=1.D0 
DO 100 J=l,9 
100 VARD ( J) =C ( 1 » J) 

DO 200 1=2.5 
T1=TTRANS*T1 
DO 200 J=1 . 9 

200 VARD (J) = VARD (J)+T1*C(I.J) 

DAY=AINT (DAY2) 

FDAY=UAY2-DAY 

THETG=THETGO ( IY1)+THD0T1*DAY+THD0T2*FDAY+EQ 
CALL EGRAV 

I F ( .NOT . ISwT(22) ) GO TO 8 
C COMPUTE LUNAR GRAVITY EFFECTS 

DO 9 1=1.4 
9 XYZ ( I ) =VARD ( I ) 

LS=1 

CALL SLGRAV 

8 IF ( .NOT • ISWT (23) ) 60 To 16 
C COMPUTE SOLAR GRAVITY EFFECTS 

DO 14 1=5.8 
IN=I-4 

14 XYZ ( IN) =VARD ( I ) 

LS=2 

CALL SLGRAV 

16 IF ( .NOT • ISWT (24) ) GO TO 18 
DO 7 1=5.7 
lN=I-4 

7 XYZ(IN)=VARD(I) 

CALL DRAG 

18 IF ( .NOT . ISWT (25) ) GO TO 20 

C COMPUTE EFFECTS DUE TO SOLAR RADIATION 
DO 19 1=5.8 
IN=I-4 

19 XYZ ( IN) = VARD ( I ) 

CALL SHADOW (IN 0 ) 

IF ( IND .EQ.O) GO TO 20 
CALL SOLRAD 
20 DO 6 1=1.3 
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6 XDD(K.I)=XXDD(I)*TT 
GO TO 25 

10 DO 15 1=1.3 

XXDD ( I ) =-X ( K . I ) /RQ*GM 
15 XDD(K.I)=XXDD(I)*TT 
25 CONTINUE 

CALL CLOCK (XR) 

FORCS=FORCS + XR - XE 
IF(R.GT.l.DO) 60 TO 30 
WRITE (3.11) 

11 FORMAT ( 1H0 » 67HRADIUS VECTOR OF SATELLITE LESS THAN EARTH RADIUS 
♦-RUN TERMINATED.) 

STOP 

30 RETURN 

end 


Q ELT HEMINT » 1 . 670830 » 1+7590 
a EOF U 

SUBROUTINE HEMINT 

N=ORDER OF INTERPOLATION POLYNOMIAL - UPPER BOUND =40. 

X=ARRAY CONTAINING N PLUS ONE EQUI-SPACEU DATA POINTS 
XD=ARRAY FOR VELOCITIES - SPECIFICATIONS AS FOR X 
H= INTERVAL BETWEEN POINTS 

T I=TIME AT WHICH POINT IS TO BE INTERPOLATED 

T= START TIME T ZERO 

FX= INTERPOLATED POSITION VALUES 

FXD= INTERPOLATED VELOCITY VALUES 

DOUBLE PRECISION FX . FXD » FXDD > TI 

DOUBLE PRECISION X * XD . XDD . XXDD r DIFF . CDEL . TT . T . TREQ » TOUT 

DOUBLE PRECISION STEPD . H. TEMP4 . OLDT » TM. TEMP3 

DOUBLE PRECISION T0L1.T0L2.TC 

DOUBLE PRECISION HX . HXB » HXB2 

DOUBLE PRECISION SMI » SJMI . SMJ » FJ * SHSQ ► COFF 

DOUBLE PRECISION ONE.TWO.THREE.VJX.ZJX.WUX 

DOUBLE PRECISION HDX ► HDD » HDBX » HDDS >HX3B» HDB8 » HDDBB 

DOUBLE PRECISION SMJSQ . SJMIQ » SMISQ » AQ » SMIQ 

DOUBLE PRECISION FN » F I . LS . NF » SI » S » SH 

DOUBLE PRECISION SMI 1 » SMI2 * FIM1 . FNI2 . A 

LOGICAL ISWT.SW 

COMMON/ I NTERP/FX (60.3) »FXD(60.3) .FXDD(b0.3) .TI.K1.M.L 
COMMON/LIMI TS/TOLl » T0L2 ► TC . ISCT r ISWT ( 30 ) > SW ( 20 ) . ORDER . LI » L2 . MODE 
COMMON/ WORKER/X (200.3) »XU(20U»3) » XDD (200.3) .XXDD (3) .DIFF (60 .3) . 

1 CDEL.TT.T. TREQ. TOUT. K.N 

C0MMON/TIMING/STEPD(90.3) .H.TEMP4.0LDT.TM.XM. YM.ZM.XS. 

♦ IENT.ITERS.INCOWL.ICH 

DIMENSION A ( 20 ) 

IF ( ISCT .GT. M) GO TO 2 
M=ISCT 

SW ( 11 ) = .TRUE • 

2 CONTINUE 

IF(M.LT.b) SW ( 11 ) = . TRUE . 

IF ( • NOT .SW ( 10 ) ) GO TO 3 

TEMP3=H 

H=-CDEL 

3 CONTINUE 
DO 1 U=1 » 3 
FX(L» J)=O.DO 
FXDIL. J)=0.DU 

1 FXDU(L> J)=0.DU 


LS=1.0D0 
NF=1 . 0D0 
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S1=0.D0 

fn=m 

S=(TI-T)/H 
SH=S *H 
ShSQ=SH**2 
DO 4 1 = 1 * M 
F 1 = 1 

SM11=S-(FI-1.0D0> 

LS=LS*SMI1 

NF=NF*FI 

4 S1=S1+1.0D0/SMI1 
NF=NF/FN 

A ( 1 ) = ( (-l.DO)**M*LS)/(NF*S) *1-1.00) 

DO 5 I=2»M 
FI = I 

5MI1=S-(FI-1.D0) 

SMI2=S-(FI-2.0D0) 

FN12=(FN-1.DU)-(FI-2.D0) 

FIM1=FI-1.0D0 

A(I)=-1.0D0*A(I-1)*(SMI2/SMI1)*(FNI2/FIM1) 

5 CONTINUE 
J=0 

F J=J 
10 1=0 

SMI=0 • QO 
SMIQ=0.D0 
SJM1=0.D0 
SJMIQ=0.U0 
LN=K1-J 
15 F 1=1 

IF (J .EQ. I) GO TO 20 
SM1 = SMI + l.UD0/(S-FI ) 

SM I Q=SM I Q+ 1 . D 0 / ( S-F I ) **2 
SJMI= SJMI + 1.000/(FJ-FI) 

SJMIQ=SJMIQ + 1.00/ (FJ-FI ) **2 
20 1=1+1 

IF (1 .LT. M) GO TO 15 
SMISQ=SM1**2 
SMJSU=SJMI**2 
SMJ=S-FJ 
AG=A ( J+i ) **3 
IF (Swill) ) GO TO 50 
C0FF=A(J+1)**2 
HX=COFF-2.DO*SMJ*SJMI*COFF 
HXB=H*SMJ*COFF 
HXB2=SMJ*C0FF/H 
25 UO 30 I 1=1 t 3 

FX(L»II)=FX(L»II) + HX*X(LN»II)+FIXB*XD(LN.II) 

30 FXD(L»II) =FXD (Lr I I ) + HX*XD ( LN. I I ) +HXB2+XD0 (LN » I I ) 
J=J+1 
F J- J 

IF ( J .LT. M) GO TO 10 

IF ( . NOT • SM ( 10 ) ) RETURN 

H=TEMP3 

Sift ( 10 ) =. FALSE • 

RETURN 

50 COFF = 3. DO * AQ/H 

ONE= 3. DO * SMJSQ + SJMlQ 
TWO= SM J -3. DO * SMJ**2 * SJMI 
THRE£= 3. DO * SMISQ - SMIQ 
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VJX= 1.00 + 3.00/2.00 * SMJ**2 * ONE - 3.00 * SMJ * SJMI 
ZJX= H**2/2.D0 * SMJ**2 
WJX= H*TWO 

HX=(1.D0 + 3.00/2.00*SMJ**2 * (3.00 * SMJSQ ♦ SJMIQ) 

* - 3. DO * SMJ*SJMI)*AQ 

HDX= (VJX* SMI + SMJ * ONE - SJMI)* COFF 

HDD= (ONE + 6. DO * SMI * SMJ * ONE - 6.00 * SMI*SJMI 

* + VJX * THREE) * COFF/H 

HXB = (SMJ -3.D0*SMJ**2*SJMI) * H * AQ 

HDBX= (3. DO * TWO * SMI + 1.00 - 6.00 * SMJ * SJMI)* AQ 
HDDB = ( -2.00 * SJMI + 2. DO * SMI -12. DO* SMJ * SJMI * SMI 

* + TWO * THREE) * COFF 
HXBB = 0 • 5D0*H**2*SM J**2* AQ 

HDBB = (3. 00/2. DO * SMJ **2 * SMI + SMJ)*H*AQ 

HDDBB= (1.00 +6.00 * SMJ *SMI +3.00/2.00 * SMJ**2 * THREE )*AQ 

DO 60 II=lr3 

FX(L»II)=FX(L»II) + HX * X(LNf II ) + HXB*XD(LN» II)+ HXBB*X0D(LN» i) 

1 /H**2 

FXO(L»II)=FXO(L*II)+HDX * X ( LN r 1 1 ) + HOBX * XD(LN*II) 

* + H08B * XDD ( LN * 1 1 ) /H**2 
FXDO(L»II)=FXOO(L»II) + HOD* X(LN*II) ♦ HDDB * XD(LN#II) 

* + HODBB * XDD(LN»II)/H**2 
60 CONTINUE 

J=J+1 

FJ=J 

IF ( J .LT. M) SO TO 10 

H=TEMP3 

SW(10)=. FALSE. 

SW(11>=. FALSE. 

RETURN 

END 

Q ELT INPUT » 1*671013* 34113 
S EOF 0 

BLOCK OATA 
C 

DOUBLE PRECISION DRAO * DTWOPI * ORSEC * THETGO ► THD0T1 * TH00T2 » GM * AE ► CS 
DOUBLE PRECISION DMIN.TC1*R.RSQ*RQ»THETG«ETIME 

DOUBLE PREC I SI ON ALT * DEN * HM * WE * ESQ *B*RH01* RH02 * RH03 * RH04 *CN* CDP * 

* C APR * EMASS * X YZ * AREA * SATMAS * CD 
DOUBLE PRECISION PSUN'CSUBR'SIGMA'F 

C 

COMMON/CONST/DRAD * DTWOP I * ORSEC * RAD ► RSEC 
COMMON/CONST 1/THETGO ( 10 ) . TH00T1 * TH00T2*0MIN»ETIME» IYBEG 
COMMON/CONST2/GM » AE * R * RSQ ► RQ * THETG ► TC 1 
COMMON/CONST3/EMASS ( 2 ) ,XYZ(4) *LS 
C0MM0N/C0NST4/PSUN ► CSUBR * SIGMA ► F 

COMMON/PERTS/ALT ( 41 ) * DEN (41*2) *HM(40*2) *WE*ESQ»B*CN*CDP*RH01»RH01 

* RH03 » RH04 * C APR * AREA * SATMAS * CD ► MO 
COMMON/CSUN/ASUN* PMOON* EMOON* TASUN ( 10 ) 

COMMON/FMODEL/CS ( 20 * 23) * INDEX1 * INDEX3 

C 

C CONVERSION FROM DEGREES TO RADIANS 
C TWO PI IN RADIANS 

DATA DRAD/.0174532925l9943296D0/*DTWOPI/6.2831853071795864D0/* 

C CONVERSION FROM SECONDS OF ARC TO RADIANS - DOUBLE AND SINGLE PREC. 

C CONVERSION FROM SECONDS OF ARC TO RADIANS 

* ORSEC/. 484813681109536D-5/ * RAD/ .0174532925/ * RSEC/.484813681E-5/ 

RIGHT ASCENSION OF GREENWICH AT JAN 0.0 FOR 1960-1969 IN DEGREES 
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DATA THETGO /98.6740065D0, 99.420934700, 99. 1822166D0 , 

. 98 * 9434986D0 * 98.7047825D0 » 99.4517117D0* 

. 99 .212993600, 98.974276500 , 98.7355595D0, 

• 99. 482489600/ f 

C MEAN ADVANCE IN RT ASC OF GREENWICH PER MEAN SOLAR OAY IN DEGREES 
C MINUTES IN 5 DAYS 

* THDOTl/.9856473354D0/,DMlN/7200.D0/ 

C GRAVITATIONAL CONSTANT TIMES MASS OF EARTH IN CANONICAL UNITS 
C SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID IN VANGUARD UNITS OF LENGTH 
DATA GM/1.D0/»AE/1.D0/ 

RATIO OF MASS OF SUN TO MASS OF EARTH 

DATA EMASSI2) /332948.55D0/ 

RATIO OF MASS OF MOON TO MASS OF EARTH 

DATA EMASS(i)/. 01230012300/ 

COMMON BLOCK /CSUN/ CONTAINS CONSTANTS NEEDED TO COMPUTE SOLAR AND 
LUNAR EPHEMERIDES 

SEMI-MAJOR AXIS OF EARTH’S ORBIT AROUND SUN IN METERS 
DATA ASUN /.1496E12/, 

SEMI LATUS RECTUM OF MOON'S ORBIT IN METERS 
. PMOON /.3847509Q2E9/, 

ECCENTRICITY OF MOON’S ORBIT 
. EMOON /. 054900489/ » 

TRUE ANOMALY OF SUN AT JAN 0.0 FOR 1960-1969 IN DEGREES 

. TASUN /3. 5727587. 2.8430634. 3.0989676. 3.3548717. 

. 3.6107759. 2.8810806. 3.1369848. 3.3928890. 

. 3.6487942, 2.9190979/ 

COMMOS BLOCK C0SST3 - CSUBR=RADIATIOS COSSTAST 

PSUN=SOLAR RADIATION PRESSURE IN DYNES/CM** 
DATA PSUN » CSUBR/4 . 5D-05 , 2.000/ 


C 


COMMON BLOCK PERTS PARAMETERS FOR ATMOSPHERIC DENSITY AND DRAG. 
ALTITUDE IN KM 

DATA (ALT(I) ,I=1,41)/O.DO» .2D+02, .40+02, .60+02, .80+02, .1D+03, 

* .12D+03. .14D+03, .16D+03, .18D+03, .2D+03. 


* .220+03. .240+03, .260+03, .28D+03, .30+03, 

* .320+03, .340+03, .36D+03. .380+03, .40+03, 

* .42D+03, .440+03, .46D+03, .48D+03, .5D+03, 

* .520+03, .54D+03, .56D+03, .580+03, .6D+03, 

* .640+03, .680+03, .72D+03, .76D+03, .8D+03, 

* .840+03, .88D+03, .92D+03, .960+03, . 1D+04/ 
ATMOSPHERIC DENSITY IN SUNLIGHT - GMS/KM**3 

DATA (DEN(I.l) ,1=1 ,41)/. 12250+13, .88910+11, .399570+10, .305920+09, 

* .19990+08, .49740+06, .2490+05, .3961D+04, 


* .12550+04, .54420+03, .27990+03, . 1596D+03, 

* .97280+02, .62190+02, .41190+02, .2805D+02, 
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c 


c 


c 


* .1953D+02*. 13850+02* .9976D+01 * .7282D+01 * 

* .53760+01* .40330+01* .3034D+01* .22990+01* 

* .1754D+01* .13490+01* .10380+01* .80430+00* 

* .6260+00* .4892D+00* .38380+00* .2388D+00* 

* .15070+00* .96580-01* .62890-01* .417D-01* 

* . 2910-01*. 20750-01 * . 15D-01 * . 111D-01* 

* .840-02/ 

ATMOSPHERIC DENSITY IN EARTH’S SHADOW 

DATA (DEN(I*2) * 1=1*41) /. 1225D+13* .8891D+11* .399570+10 . .30592D+09, 


* 

* 

* 

* 

* 

* 

* 

* 

* 


• 1999D+08* ,49740+06* .2490+05* .4104D+04* 

. 12820+04*. 5109D+03* .23270+03* . 11590+03* 
.61740+02* .34670+02* ,20310+02* .1231D+02* 
. 76750+01*. 4894D+01* .31780+01*. 2094D+01* 
. 13970+01 * . 96250+00 * . 65520+00 * . 44970+00 * 
.31090+00* .21660+00* .152D+00* .10760+00* 

. 76850-0 1 * . 5545D-0 1 . . 40470-0 1 * . 22440-01 * 
.13250-01* .83970-02* .57130-02* .4143D-02* 
.3150-02* .2520-02* .208D-02* . 1765D-02* 


* .1510-02/ 

0RAG CONSTANTS 

DATA WE *CN * CAPR * CD/ . 588336903074954198Q-0 1 * 298 . 2500 * 

* 6378 . 166D0 * 2 • 300/ 

AREA AND MASS OF SATELLITE 
DATA AREA * SATMAS/ . 914838700+04 * . 635029350+05/ 


DIFFERENTIAL CORRECTION PARAMETERS AND ANGLE OF ATM. BULDGE 
DATA RH01 * RH02 * RH03 * RH04/ 3*U .00 * 30 .DO/ 

DATA MD/6/ 

COMMON BLOCK /FMODEL/ CONTAINS THE GRAVITATIONAL COEFFICIENTS AND 
THEIR INDEXES 


INDEXES FOR COEFFICIENTS 
OATA IN0EX1*IN0EX3 /16* 16/ 


•C ( 2 * 0 ) THRU C (20*0) 

OATA (CS(I*1) *1=2*20) /-1082. 6450-6*2. 5460-6*1. 6490-6*. 210D-6* 

-.6460-6* .3330-6* .2700-6* .5300-7* .540D-7* 
-.3020-6* .3570-6* .1140-6* .1790-6*6*0.00/ 

•C ( 2 * 1 ) THRU C (20* 1) 

DATA (CS(I*2) *1=2*20) /0. 00*2. 0910-6*-. 5430-6*-. 677D-7*-.037D-6* 

. 1440-6*-. 052D-6* .0760-6* .6490-7* -.3130-7, 
, -.9230-7 *0.00 *-.7880-8 *6*0. 00/ 

■C(2*2) THRU C(20*2> 

OATA (CS(I*3) *1=2.20) /l. 5360-6* .2510-6* .7380-7* .1020-6* .8580-8* 

, .3630-7* .21350-8*-. 2770-9*-. 6240-8*0. 00* 

, -.4900-8*8*0.00/ 


C(3*3) THRU C(20*3) 

DATA (CS(I*4) .1=3*20) / .7820-7* . 5090-7*-. 1720-7*-. 1120-8* . 3520-8. 

-.3740-9* 0.00* -.3790-9. 10*0.00/ 


C(4»4) THRU C (20 *4) 

DATA (CS(I*5) *1=4*20) /-. 1120-8*-. 2060-8*-. 1670-9*-. 323D-9* 
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-. 2770-9 * 0. DO 4360-10 * 10*0. DO/ 

-C(5*5> THRU C (20 * 5 ) 

DATA (CS ( I * 6) *1=5*20) /.384D-9*-.253D“9* .2690-10 *-.9590-11 *12*0./ 
*/ 

— C (6*6) THRU C(20*b) 

DATA (CS( I *7) * 1=6*20) /-. 932D-11 *-, 1450-10 * -.4750-12 * 12*0.00/ 
-C(7*7) THRU C(20*7) 

DATA (CS ( I * 8) *1=7*20) /. 1020-11*-. 4440-13.12*0. 00/ 

-C(8*8) THRU C (20 ► 8) 

DATA (CS(I*9) .1=8*20) /-. 316D-12 * 12*0 . DO/ 

-C (9*9) THRU C(20*9) 

DATA (CS(I*10) *1=9*20) /6*0.D0*-.24lD-18*5*0.U0/ 

-C(10*10) THRU C (20 * 10 ) 

DATA (CS(I*11) .1=10*20) /11*O.QO/ 

-C(ll*ll) THRU C(20*ll) 

DATA ( CS ( I .12) * 1=11*20) /3*0. DO * .9470-21 >6*0. DO/ 

-C ( 12* 12) THRU C (20 * 12) 

DATA (CS( 1*13) *1=12.20) /-. 27830-18. -.1100-18* .504D-19*-. 1140-19, 

5*0. DO/ 

-C(13*13) THRU C ( 20 * 13) 

DATA (CS(I*14) *1=13*20) /-.216D-19*0. 00*-. 1170-20*5*0. 00/ 

“C ( 14 * 14 ) THRU C (20 * 14) 

DATA (CS(I*15) *1=14*20) /-. 1931D-21 * . 114D-22*5*0.D0/ 

-C ( 15* 15) THRU C (28* 15) 

DATA (CS(I*16) *1=15*20) /6*O.DO/ 

-C ( 16* 16) THRU C (20 * 16) 

OATA (CS(I*17) *1=16*20) /5*O.DO/ 

-C ( 17* 17) THRU C (20 * 17) 

DATA (CS(I*18) .1=17*20) /4*O.DO/ 

-C ( 18*18) THRU C ( 20 * 18) 

DATA (CS(I*19) * 1=18. 20)/3*0.00/ 
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•C(19# 19) THRU C (20* 19) 

OATA (CS(I»20) #1=19#20) /2*0.D0/ 

■C (20 # 20 ) 

DATA CS ( 20 # 21 ) /O.DO/ 

■S<20#0) THRU S(2»0) 

OATA (CS(I#23> #1=1.19) / 19*0.00/ 

•S (20 » 1) THRU S(2»l) 

OATA (CS(I#22) #I=1#19) /fa*0D0# . 2800-8.0. 00#-. 0400-6# .0880-7# 

-.7790-7# .7840-8# .045D-6# .1140-6#-. 212D-7, 
-. 8820-7. -.445D-6# .287D-6#0.D0/ 


•S(20#2) THRU S ( 2 #2 ) 

DATA (CS(I#21) #1=1 #19) /8*0. 00#-. 023D“8#0. 00#-. 2500-8# .242D-8# 

.3200-8# .162D-7#-. 4550-7#-. 375D-7, .1480-6, 
- . 1 840-6 # - . 8720-6/ 

■S120.3) THRU S(3#3) 

DATA (CS ( I #20) # 1=1 #18) /10*0 .00# . 175D“9# 0 .00 # .4040-10 # .254D-9# 

.6430-9# .2310-9#-. 1140-7# .2260-6/ 

■S (20 # 4) THRU S(4#4) 

OATA (CS ( I » 19) » 1=1 » 17) /10*0.D0 #-. 6540-10 # 0.D0#-. 1570-10 # -.2170- 9 . 

-.1960-8# , 498D-9# .4860-8/ 


•S (20 * 5) THRU S(5#5) 

DATA (CS ( I # 18) ► 1=1 # 16) /12*0 .DO *. 2140-10 #. 1910-10 #-. 3700-9# 

-.1460-8/ 


S(20#6) THRU S ( 6# fa ) 

OATA (CS ( 1 # 17) » 1=1 ► 15) /12*0. DO# .8880-11# .4730-11#-. 3610-10/ 
S(20>7) THRU S(7»7) 

DATA (CS(I#lb) #I=1#14) /12*0. 00# .1580-12# .0180-10/ 

S(20#8) THRU S ( 8 # 8 ) 

OATA (CS ( I « 15) # 1=1 # 13) /12*0. 00# .1300-12/ 

S(20#9) THRU S(9#9) 

DATA (CS(I»14) #1=1 #12) /5*0. 00 #-.4830-18# 6*0. 00/ 

S (20 # 10) THRU S(10»10) 

OATA (CS(I»13) #1=1*11) /11*0.D0/ 

S(20# 11) THRU S(11#1D 
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DATA (CS ( I * 12) *1=1*10) /6*0 . DO *-. 473D-21 * 3*0 . DO/ 

>5(20*12) THRU 5(12*12) 

DATA (CS(I»11) *1=1*9) /5*0.D0*.l068D-19*-.150D-19*.933D-19* 

.718D-20/ 

■S(20*13) THRU S( 13* 13) 

DATA (CS(I*10) *1=1*8) /5*0. DO*-. 928D-21*0, DO* .2820-19/ 
■S(20* 14) THRU S ( 14 * 14 ) 

DATA (CS(I*9) *1=1*7) /5*0.D0*-.558D-22*-. 41380-22/ 

•S(20*15) THRU S(15*1S) 

DATA (CS( I * 8) * 1=1 .6) /6*0.D0/ 

S (20 * 16) THRU S ( 16 * 16 ) 

DATA (CS ( I r 7) *1=1*5) /5*0.D0/ 

5(20*17) THRU S ( 17* 17) 

DATA (CS(I*6) *1=1*4) /4*0 .00/ 

S(28* 18) THRU 5(16*18) 

DATA (CS(I*5) *1=1*3) /3+0.D0/ 

S (20 ► 19) THRU S ( 19 * 19) 

DATA (CS(I*4) *1=1*2) /2*0.D0/ 

5(20*20) 


DATA CS(1*3) /0.D0/ 

C 

END 

a ELT OUTPUT * 1 *670830 * 47591 

a eof a 

SUBROUTINE OUTPUT 

DOUBLE PRECISION TOUT * TX * TXD *EV * EVD* X * XD ► XDD * XXDD* DIFF*CDEL* TT * T 
10UT *T0L1*T0L2*TC*BEGTIM*ENDTIM*BEGT2*DSQRT*TREQ 
DOUBLE PRECISION FX*FXD*FXDD*TI*TIME 
DOUBLE PRECISION DRAG*DN 
DOUBLE PRECISION RSAVE 
LOGICAL ISWT * SW 
INTEGER ORDER 
DIMENSION TX ( 3) * TXD ( 3) 

COMMON/ WORKER/X (200*3) *XD(200*3) * XDD (200 *3) * XXDD (3) *DIFF(60*3) * 

1 CDEL*TT*T*TREQ*TOUT*K*N 

COMMON/ I NTERP/FX (60*3) *FXD(60*3) *FXDD(60*3) *TI*K1*M*M1 
COMMON/LIMI TS/TOL1 * T0L2 * TC * ISCT * ISWT ( 30 ) * SW ( 20 ) * ORDER * LI * L2 * MODE 
COMMON/OPT/BEGT IM * ENDT IM , 8EGT2 
COMMON/DDRAG/DR AG ( 3 ) 

COMMON/BRIEF/RSAVE ( 10 ) 

101 FORMAT(1HO*21HNORM OF ERROR VECTORS * 8X2HE=D24. 16* 5X5HED0T=D24. 16 ) 



67 


107 FORMAT ( 1HQ . 15X1HT27X1HX30X1HY28X1HZ ) 

108 FORMAT (4(5X»D24.16) ) 

109 FORMAT ( 14X5HR X V . 15X . 3 (024. 16. 5X) ) 

110 FORMAT ( 1H0 . 14X . 20HNORM OF DRAG VECT0R=024 . 16) 

TOUT=T/TC 

IF(.NOT. ISWT (7) ) GO TO 57 

1FIT0UT.LT.BEGTIM.0k. (TOUT. GT.ENOTIM. AND. TOUT. LT.BEGT2) ) RETURN 
57 WRITE ( 3. 107) 

IF (ISWT (6) ) 60 TO 20 

WRITE (3. 108) TREQ . (FX (Ml » I ) .1=1.3) 

CON=DSQRT( (FX(M1.2)*FXD(M1.3)-FX(M1.3)*FXD(M1»2> )**2 

* + (FX (Ml . 3) *FXD(M1 . 1 ) -FX (Ml » 1 ) *FXD (Ml » 3) >**2 

* +(FX(M1.1)*FXD(M1.2)-FX(M1.2)*FXD(M1.1) )**2) 

WRITE (3. 108) CON. (FXD (Ml . I ) . 1=1 . 3) 

IF (ISWT ( 12) ) GO TO 26 

ON=O.DO 

DO 23 1=1.3 

23 DN=AMAX1 (DABS (DR AG ( I ) ) »DN) 

IF(DN.LE.O.DO) GO TO 24 
WRITE( 3. 110) DN 

24 CONTINUE 
RETURN 

26 T IME=T 
T=TI 

CALL TRANS(TX.TXD) 

T=TIME 

ev=o.do 

EVD=O.DO 
DO 22 1=1.3 

EV=AMAX1 (DABS (FX (Ml . I ) “TX ( I ) ) .EV) 

22 EVD=AMAXl(DABS(FXD(Ml.I)-TXD(I) ) .EVD) 

RSAVE ( 3) =£V 
WRITE ( 3. 101 ) EV .EVD 
RETURN 
20 CONTINUE 

WRITE (3. 108) TOUT. ( X (K . I ) . 1=1 » 3) 

WR1TE(3.109> (XD(K.I) .1=1.3) 

CON=OSQRT( (X(K.2)*XU(K.3)-X(K.3)*XD(K.2) ) **2 + ( X (K. . 3) *XD ( K . 1 ) -X (K 
1)*XD(K»3) )**2+(X(K.l)*XD(K»2)-X(K.2)*XD(K.l) )**2) 

WRITE (3.108)CON .XXDD 
IF (.NOT. ISWT ( 12) ) RETURN 
CALL TRANS(TX.TXD) 

25 EV=O.DO 
£VD=0 »D0 

DO 32 1=1.3 

EV=AMAX1 (DABS (X (K. I ) “TX ( 1 ) ) .EV) 

32 EVD=AMAX1 (DABS (XU(K.I)-TxD(I) ) .EVD) 

KSAVE(3)=EV 

WkITE(3.101)EV.EVD 

return 

END 

N ELT PCCOFF. 1.670830. 47592 
M EOF a 

SUBROUTINE PCCOFF (SIGMA.GAMMA. SIGMAS. gammas. M) 

C CALCULATES PREDICTOR CORRECTOR COEFFICIENT M GIVEN M-l COEFFICIEVT 

DOUBLE PRECISION SIGMA. GAMMA. SIGMAS. GAMMAS. SI. S2. S3. S4.H. Xl .X2.FM, 
1FI 

DIMENSION SIGMA (63 ) .GAMMA (63) »SIGMAS(63) .GAMMAS(63) »H(63) 
lF(M-2) 1.2.3 
1 SIGMA(M)=1.0DU 
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SIGMAS (M)=1.0D0 
GAMMA(M)= 1.000 
GAMMAS (M)=1.0D0 
H(M)= 1.0D0 
GO TO 10 

2 SIGMA(M)=0.000 
SIGMAS (M ) - -1.000 
GAMMA (M) =0 • 500 
GAMMAS (M)=-0.5D0 
H(M)= 1.5D0 

GO TO 10 

3 S1=0.000 
S2=0.0D0 
S3=0.0D0 
S4=0.0D0 
FM=M 

H(M)=H(M-1)+1.000/FM 

J=M-1 

DO 4 I=1»J 
FI-I 

Xl=l.ODO/(FI+l.OD0) 

X2=2.0D0/(FI+2.OD0) 

K=M-I 

L=I+1 

S1=S1+X1*GAMMA(K) 

S2=S2+X1*GAMMAS(K) 

S3=S3+X2*H(L)*SIGMA(K) 

4 S4=S4+X2*H(L)*SIGMAS(K) 

GAMMA ( M ) =1 . 0D0-S1 
GAMMAS(M)=-S2 

SIGMA (M)=1.0D0-S3 
SIGMAS (M) = -S4 
10 RETURN 
END 

» ELT RESUME » 1 f 670830 > 47593 
(i EOF M 

SUBROUTINE RESUME 
DOUBLE PRECISION T0L1 f T0L2 r TC 
DOUBLE PRECISION RSAVEfTEMP 
LOGICAL ISWT t SW 

DOUBLE PRECISION STEPD*TEMP3>TEMP4rOLDT f TM 
DOUBLE PRECISION XfXQfXDDfXXDDfDIFFfCDELfTT fT fTREQfTOUT 
COMMON/WORKER/X (200 r 3) fXD(200f3) fXDD(200f3) fXXDD(3) f0IFF(60f3) f 
1 CUEL f TTf T f TREQ fTOUT f K f N 

COMMON/LIMI TS/T0L1 f T0L2 f TC f ISCT f ISWT ( 30 ) > SW ( 20 ) f ORDER f LI f L2 f MODE 
COMMON/TIMING/STEPD (90.3) • TEMP3f TEMP4 fOLUT f TM f XM f YMf ZM f XSf 
* IENTfITERSfINCOWLfICH 

COMMON/ TIMES/ TAB f COWL fFORCS 
COMMON/BRIEF/RSAVE ( 10 ) 

C INSERT CODING FOR RESUME 

IF (.NOT. SW ( 13) ) GO TO 1 
IF (M0DE.EQ.3 .OR. MODE. EG. 4) GO TO 58 
TEMP3=TEMP4 
1 J=1 

92 IF(DA6S(TEMP3-STEPD(Jf 1) ) .GT. 1.0-13) GO TO 90 
STEPD ( J r 2 ) = STEPD(Jf2) + DBLE(ZM) 

STEPD ( J r 3) = STEPDt J»3)+ TM 
IF ( SW l 13) ) GO TO 58 
RETURN 
90 J=J+1 
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IF(J.LE.IENT) GO TO 92 
IF(IENT.GT.89) GO TO 190 
1ENT= IENT + 1 
STEPD ( IENT » 1 ) =TEMP3 
STEPO ( IENT . 2 ) =DBLE ( ZM ) 

ST£PD( IENT .3) = TM 
IF(.NOT. SW ( 13) ) RETURN 
58 ZM=YM-XS 
XS=YM 

TAB=TAB/60.D0 

COWL=COWL/60.DO 

FORCS=FORCS/60.DO 

ZM=ZM/60. 

WRITE(3.360) TAB.COwL.FORCS.ZM 

360 FORMAT (1H1. //////» 47X. 27HTIMING BREAKDOWN BY ROUTINE //. 

1 17X. 5HTABLE * 24X» SHCSTEP. 24X. 4HFRCS » 25X. 

2 9HT0TAL RUN /(4(5X»E24.8> ) ) 

WRITE(3»362) (RSAVE ( I ) * 1 = 1 ► 3) 

362 FORMAT (1H0 // 32X. 13HINITIAL ORDER 10X. 12HINITIAL STEP 10X. 

* 11HFINAL ERROR / 33X. 3(D12.5» 10X) ) 

AITER = FLOAT(ITERS)/FLOAT(INCOWL) 

WRITE (3.365) AITER 

365 FORMAT (1H0 //// 41X. 37HAVERAGE NUMBER OF ITERATIONS IN CSTEP// 

1 48X. E24.8) 

WRITE(3.366> INCOWL 

366 FORMAT (1H0//// 50X. 28HT0TAL NUMBER OF COWELL STEPS//. 57X. 112) 
IF(MOOE.EQ.l .OR. M0DE.EQ.2 .AND. IENT.GT.l) GO TO 76 
STEPD(1»1)=CDEL/TC 

STEPD (1.2) =DBLE ( ZM ) 

STEP0(1.3)=T/TC 

WRITE(3»370) (STEPD ( 1 ► I ) .1=1.3) 

RETURN 

76 IT=IENT-1 

DO 80 11=1. IENT 
STEP0(H.1)=STEPD( Il.D/TC 
STEPU(I1.2)=STEPD(I1.2)/60.00 
80 STEPD ( II . 3)=STEPD ( II » 3) /TC 

78 DO 77 1=1. IT 
12=1 

DO 77 LR=I2. IT 

IF (STEPD (LR+1 » 1 ) .GT . STEPD ( I » 1 ) ) GO TO 77 

DO 79 J=1 » 3 

TEMP= STEPD (LR+l.J) 

STEPD (LR+1 .J)=STEPD( I » J) 

79 STEPD ( I ► J) =TEMP 

77 CONTINUE 

WRITE (3.370) ( ( STEPD ( I ► J ) . J=1 . 3 ) . 1=1 . IENT ) 

370 FORMAT (//////. 22X. 9HSTEP SIZE 27 X. 8HRUN TIME. 21X. 

* 10HORBIT TIME (//. 20X . D24 . 16 . 5X . D24 . 16. 5X . 024 . 16) ) 

RETURN 

C ERROR EXIT 

190 IF ( .NOT • SW( 13) ) RETURN 
WRITE (3.290) 

290 FORMAT (49H0STEP CHANGES GREATER THAN 90. RESUME INCOMPLETE.) 
RETURN 
END 

0 ELT RK.l. 670830. 47596 
0 EOF Q 

SUBROUTINE RK 

C SHANKS R-K ROUTINE ORDER 8 
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DOUBLE PRECISION X * XD . XDD * H * T * T I ME * F * RK 1 » RK2 * RK3 * RK4 r RK5 * RK6 * RK7.X 
11*XD1*RK8*RK9*RK10*XXDD*DIFF*CDEL*TT *T0L1*T0L2*TC*TREQ*T0UT 
DOUBLE PRECISION H1.H2 
DOUBLE PRECISION CSTEPT 
INTEGER ORDER 
LOGICAL ISwT *SW 

DIMENSION F(6) *RKl(6) *RK2(6) *RK3(6) *RK4(6) *RK5(6 

1) »RK6(6) * RK7(6) *RK8<6) *RK9(6) *RK10(b) rXl(3) *XD1(3) 

COMMON/ WORKER/X (200*3) rXD<200f 3) »XDD(200»3) »XXDD(3) »DIFF(60»3) » 

1 CDEL * TT»Tr TREQ* TOUT » K » N 

COMMON/RKT/CSTEPT 

COMMON/LI M I TS/TOL1 » T0L2»TC ► ISCT r ISWT ( 30 ) » SW (20 ) r ORDER r LI >L2r MODE 

COMMON/RKST/HI rH2 

IF (SW(1) ) GO TO 110 

H=-H2 

LL=3 

DO 73 1=1 » 3 
X(Kt I)=XCK+1* I) 

73 XD (K * I ) =XD (K+l 1 1 ) 

82 IF (CSTEPT «LE. (T+H) ) GO TO 120 
GO TO 61 
110 LL=1 
H=H1 

IF (SW(12) )H=H2 
DO 74 1=1*3 
X(K*I)=X(K-1*I) 

74 XD(K*I)=XD(K-1»I) 

130 IF (CSTEPT. GT. (T+H) ) GO TO 120 
61 LL=2 

H=CSTEPT-T 
120 TIME=T 

DO 1 1=1*3 
XI ( I ) =X (K * I ) 

1 XD1 ( I ) =XD (K* I ) 

DO 50 L=1 * 10 
DO 2 1=1*3 

F ( I ) =XD (K* I ) 

2 F ( 1+3) =XXDD < 1 ) 

GO TO (10*20*30*40*52*60*70*80*90*100) *L 
10 DO 3 1=1*6 

3 RK1 ( I ) =H*F ( I ) 

DO 4 1=1*3 

X (K* I ) =X1 ( I ) +RK1 ( I ) *4. 0D0/27 . ODO 

4 XD ( K * I ) =XD1 ( I ) +RK1 ( 1+3) . 000/27, 0D0 
T=TIME+H*4.0D0/27.0D0 

GO TO 49 
20 DO 5 1=1*6 

5 RK2 ( I ) =H*F ( I ) 

DO 6 1=1*3 

X(K*I)=X1(I)+(RK1(I)+3.0D0*RK2(I) J/18.0D0 

6 XD (K * I ) =XD1 (I) + (RK 1(1+3) +3. ODO+RK2 (1+3) )/18.0D0 
T=TIME+2. 0D0*H/9, 000 

GO TO 49 
30 DO 7 1=1*6 

7 RK3( I )=H*F ( I ) 

DO 8 1=1*3 

X(K*I)=X1(I)+(RK1(I)+3.0D0*RK3(I) )/12.0D0 
B XD(K * I )=XD1 (I) + (RK 1(1+3) +3. 0D0*RK3( 1+3) ) / 12.000 
T=TIME+H/3.0d0 
GO TO 49 
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40 DO 9 1=1*6 
9 RK4 ( I ) =H*F ( 1 ) 

DO 12 1=1*3 

X(K*I )=X1(I) + (RK.1(I)+3.0D0*RK4(I) 1/8.000 

12 XD(K* I)=XD1(I)+(RK 1(1+3 ) +3. 0D0+RK4 ( 1+3) )/8.0D0 
T =TIME+H/ 2.000 

60 TO 49 
52 DO 13 1=1*6 

13 RK5 ( I ) =H*F ( I ) 

DO 14 1=1*3 

X(K* I >=X1 m + (13.0D0*RKl( I)-27.0D0*KK3( I >+42.0D0*RK4(I)+8.0D0*RK5‘( 
11))/54.0DO 

14 XD(K* I)=XDl(l) + (13.l)D0*RKl(I + 3)-27.0D0*RK3(I+3)+42.0D0*RK4(I+3)+S'. 
lDU*RK5(I+3))/54.0U0 

T=TIME+2.0DO*H/3.0DU 
GO TO 49 
60 DO lb 1=1*6 
lb RK6 ( I ) =H*F ( 1 ) 

DO lb 1=1*3 

X(K*I)=Xl<I)+(389.0D0*RKl(I)-b4.0D0*RK3(I>+966.D0*RK4(I)-824.D0*fiK 
lb ( I ) +243. ODO*RKb ( 1 ) ) /4320 . ODO 

16 XD(K*I)=XDl(I> + (389.0D0*RKl(l + 3)-54.D0*RK3(l+3)+966.D0*RK4(I+3)-:?-- 
14 . ODO*RKb (1+5) +243* UuO*RK6(I+3) )/4320.0D0 

T=TIME+H/6, ODO 
GO TO 49 
70 DO 17 1=1*6 

17 RK7 ( I ) =H*F ( I ) 

DO 21 1=1*3 

X(K*l)=Xl(I)+(-231.ODO*RKl(I>+8l.ODO*RK3(I)-1164.DO*RK4(I)+b56.D0+ 
1RK5( I ) -122, OuO*HKb (I)+800.0D0*RK7(I) )/20.0D0 

21 XD<K*I)=XDl(I>+(-231.0DO*RKl<I+3)+81.DO*RK3(I+3)-1164.DO*RK4(I+3)+ 
1656. GQ0*RKb( 1+3) -122. 0D0+RK6 ( I+3)+800.0D0*RK7( 1+3) )/20.0D0 

T=TJ ME+H 
GO TO 49 
80 DO 22 1=1*6 

22 RK8( 1 ) =H*F < I ) 

DO 23 1=1*3 

X(K* I)=Xl(I)+(-127.0DO*RKl(I)+18,ODO*RK3(I>-678.DO*RK4(I)+456.DO** 
lKb(I)-9.OD0*RK6(I)+b7b.0D0*RK7(I)+4.OD0*RK8(I) )/288.0D0 

23 XD(K.*I)=XDl(I) + (-l 27. 0D0*RKl(l+3)+18.D0*RK 3(1 + 3) -678. D0*RK4 (1+3)4^ 
lbb.0D0*RKb(I+3)-9.0D0*RK6( 1+3) +576. 0D0+RK7 ( 1+3) +4. 0D0+RK8 ( 1+3) )/ 2 ? 
28. ODO 

T=TIME+5. 0D0+H/6.0D0 
GO TO 49 
90 DO 24 1=1*6 

24 KK9( I ) =H*F ( I ) 

DO 2b 1=1*3 

X (K* I ) =X1 ( I ) + ( 1481 « 0D0+RK.1 (I)“81.0D0*RK3(I) +7104.D0+RK4 ( X ) -3376 
l+RK5(I)+72.0D0*RKb(l)-5040.0D0*RK7(l)-60.0D0*RK8(I)+720.0D0*RK9(I; 
2)/820.0D0 , 

25 XD (K*I)=XD 1(I)+(1 481, ODO+RK 1(1+3) -81, D0*RK3(I+3) +7104 .D0+RK4 ( 1+3)- 
13376. 0D0+RK51 1+3) +72 . 0D0+RK6 ( 1+3) -5040 . 0D0*RK7 ( 1+3) -60 . 0D0+RK8 ( 1+3 
2) +720 . 0D0+RK9 (1+3) )/620.0D0 

T=TIME+H 
GO TO 49 
100 DO 26 1=1*6 

26 RK10 ( 1 ) =H*F ( I ) 

GO TO 51 

49 CALL FRCS 

50 CONTINUE 
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51 DO 18 I=l»3 

X(K»I)=X1(I) + (41.0U0* (RKl ( I ) +RK10 ( I ) ) +27. 0D0* (RK4 (I ) +RK6( I ) ) +27C. 
10D0*RK5 ( I ) +216. 0D0* (RK7 ( I ) +RK9( I ) ) )/840.0D0 
18 XD(K»I)=XD1(I) +(41.0D0*(RKl(I+3)+RK10(I+3) ) +27. 0D0* (RK4 ( 1+3) +RK& ( 
11+3) ) +272. 0D0*RK5 ( 1+3) +216. 0D0* (RK7 ( 1+3) +RK9 (1+3) ) )/840.0DO 
CALL FRCS 

00 TO ( 130 *81*62) rLL 
81 RE7URN 
END 

ELT RYMDI » 1 > 670830 r 47596 
EOF 8 

NAME SUBROUTINE RYMDI 

LANGUAGE FORTRAN IV 

MACHINE UNIVAC 1107/1108 


PURPOSE TO SEPARATE PACKED SIX-DIGIT DECIMAL DATES INTO 

TWO-DIGIT YEAR ► MONTH* AND DAY 


CALLING SEQUENCE 
SYMBOL TYPE 

YMD(l) I 

Y( 1) I 

M ( 1 ) 1 

D ( 1 ) I 


CALL RYMDI (YMDr Y» M» D) 
DESCRIPTION 

INPUT - DATE TO BE SEPARATED 
OUTPUT - TWO-DIGIT YEAR 
OUTPUT - TWO-DIGIT MONTH 
OUTPUT - TWO-DIGIT DAY 


ROUTINES REQUIRED NONE 
TAPES REQUIRED NONE 
CARDS REQUIRED NONE 


SUBROUTINE RYMDI (YMD» Y »M»D) 

PURPOSE - TO UNPACK YYMMDD TO YY - MM - DD 
INPUT YMD = YEAR MONTH DAY PACKED 
OUTPUT Y = YEAR 

M = MONTH 

D = DAY 

INTEGER YMD» Y »M»D 
C 

Y=YMD/10000 

I=YMD/100 

M=I-Y*100 

D=YMD-I*100 

RETURN 

END 

8 ELT SHADOW r lr671009r 46541 
8 EOF 8 
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SUBROUTINE SHADOW ( IND ) 

DOUBLE PRECISION X . XD . XDD . XXDD , DIFF . CDEL . TT . T » TREQ * TOUT 
DOUBLE PRECISION EMASS r XYZ. COSPSI . PROJ . REARTH . DSQRT 
DOUBLE PRECISION PSUN . CSUBR » SIGMA . F 
DOUBLE PRECISION GM» AE » R . RSQ. RQ . THETG. TCI 

COMMON/ WORKER/X (200.3) > XU (200. 3) » XDD (200 *3) .XXDD (3) » DIFF (60 *3) ► 

* CDEL» TT » T » TREQ . TOUT . K . N 
COMM ON/C 0NST3/ EM ASS ( 2 ) .XYZ(4) .LS 
C0MM0N/C0NST4/PSUN * CSUBR . SIGMA » F 
C0MM0N/C0NST2/GM . AE » R » RSQ . RQ . THETG . TC 1 

C COMPUTE DOT PRODUCT - GET ANGLE BETWEEN SUN AND SATELLITE 

COSPSI = ( X(K»1)/R)*XYZ(1)+(X(K.2)/R)*XYZ(2)+(X(K»3)/R)*XYZ(3) 

C DETERMINE WHETHER SATELLITE IS IN SUNLIGHT 
IF (COSPSI .LT. 0.D0) GO TO 10 
IND=1 
RETURN 

10 PROJ=DSQRT( (X(K.2)*XYZ(3)-X(K.3)*XYZ(2) )**2 

* +(X(K»3)*XYZ(1)-X(K.1)*XYZ(3))**2 

* +(X(K»1)*XYZ(2)-X(K.2)*XYZ(1) )**2) 

REARTH=1 • DO-F* l (X(K.3)+XYZ(3) * DSQRT (RSQ-1 ,00 ) )**2) 

IND=1 

IF (PROJ .LT. REARTH) IND=0 

RETURN 

END 

0 ELT SLGRAV. 1.671013. 34112 
0 EOF ID 

SUBROUTINE SLGRAV 

DOUBLE PRECISION EMASS. XYZ . GM » AE » R . RSQ. THETG. TCI » RSAT . RS 
DOUBLE PRECISION X.XD* XDD .XXDD .DIFF .COEL.TT . T .TREQ .TOUT * DSQRT 
DOUBLE PRECISION DXX(3) 

COMMON/ WORKER/X (200.3) .XD(200.3) .XDD(200.3) .XXDD(3) .DIFF(60.3) . 

* CDEL.TT.T.TREQ.TOUT.K.N 
COMMON/CONST2/GM. AE.R.RSQ.RQ. THETG .TCI 
C0MM0N/C0NST3/EMASS ( 2 ) .XYZ(4) .LS 

DO 10 1=1.3 

10 XYZ(I)=XYZ(I)*XYZ(4) 

RSAT=( (X(K.l)-XYZ(l) ) **2+ ( X ( K . 2 ) -XYZ ( 2 ) ) **2+ ( X (K . 3) -XYZ ( 3) )**2> 

RSAT=( DSQRT (RSAT) ) *RSAT 

RS=XYZ(4)**3 

DXX(1)=GM*EMASS(LS)*( ( ( X (K . 1 ) -XYZ ( 1 ) ) /RSAT ) + ( XYZ ( 1 ) /RS) ) 
DXX(2)=GM*EMASS<LS)*( ( ( X (K . 2) -XYZ ( 2 ) ) /RSAT) + ( XYZ ( 2 ) /RS) ) 
DXX(3)=GM*EMASS(LS)*( ( < X (K . 3) -XYZ ( 3) ) /RSAT ) + ( XYZ ( 3) /RS) ) 

DO 3 1=1.3 

3 XXDD ( I ) =XXDD ( I ) “DXX ( I ) 

RETURN 

END 

ID ELT SOLRAD. 1.671009. 48542 
13 EOF IS 

SUBROUTINE SOLRAD 

DOUBLE PRECISION X. XD. XDD. XXDD. DIFF. CDEL. TT.T. TREQ. TOUT 
DOUBLE PRECISION EMASS. XYZ 

DOUBLE PRECISION PSUN.CSUBR. SIGMA. F.RSS. DSQRT 
DOUBLE PRECISION S0L(3) 

COMMON/WORKER/X (200 . 3) .XD1200.3) .XDD (200. 3) .XXDD (3) . DIFF (60. 3) . 

* CDEL.TT.T.TREQ.TOUT.K.N 
COMMON/CONST3/EMASS ( 2 ) .XYZ(4) .LS 
COMMON/CONST4/PSUN . CSUBR .SIGMA . F 

DO 1 1=1.3 

1 XYZ(I)=XYZ(I)*XYZ(4) 

RSS=DSQRT ( (XYZ(l)-X(K.l) )**2+ (XYZ(2)-X(K.2) )**2 
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* + (XYK3)-X(K»3) )**2) 

00 2 1=1.3 

SQL(I)=SIGMA*< (XYZ(I)-X(K.I) )/RSS) 

2 XXDO ( I ) = XXDD ( I ) -SOL ( 1 ) 

RETURN 

END 

8 ELT SUMS .1.670830. 47597 
8 EOF 8 

SUBROUTINE sums 

double precision x.xd.xdd.xxdd.oiff.cdel.tt.t .treq.tout 

DOUBLE PRECISION S1.S2.PX.PXD.CX.CXD.CT0L.SAVE 
COMMON/WORKER/X(200.3> »XD(200»3) .XDD1200.3) .XXDD(3) .DIFF(60.3) . 

1 CDEL . TT » T . TREQ . T OUT . K . N 

COMMON/COWS/SI (3) .S2(3) »PX(63) .PXD(63) .CX(63) .CXD(63) »CTOL»SAVEd0 
1 .3) .ITER 

DO 1 1=1.3 

SI ( I >=XD (K . I ) *CDEL + CXD(2)*XDD(K» I) 

S2(I)=X(K.I)- CX ( 3) * XDD(K.I) 

DO 1 J=1.N 

S1(I)=S1(I) - CXD(U+2> * DIFF(J.I) 

1 S2 ( I ) =S2( I ) -CX( J+3) * DIFF(J.I) 

DO 2 1=1.3 

S1(I)= XDD(K.I) + SKI) 

2 S2 ( I ) = SKI) + S2 ( I ) 

RETURN 

END 

8 ELT TABLE. 1.670830. 47598 
8 EOF 8 

SUBROUTINE TABLE 

DOUBLE PRECISION X. XU . XDD » XXDD.TOUT . T . TIME. TOUT 1 .TT.D IFF .TOL1 . T0L2 
1 . CSTEPT . F J . TC . CDEL . TREQ 
DOUBLE PRECISION 6M. AE.R.RSQ.RQ.THETG.TC1 
DOUBLE PRECISION CS.UAY1. CENTER. DSTART 

integer order 

LOGICAL ISWT.SW 
COMMON/RKT/CSTEPT 

COMMON/ WORKER/X (200.3) » XU 1200.3) »XDD<200.3) .XXDD(3) »DIFF(60»3) . 

1 CDEL. TT.T. TREQ. TOUT. K.N 

COMMON/LIMI TS/T0L1.T0L2.TC. ISCT. ISWT(30) »SW(20) .ORDER. LI. L2. MODE 
DOUBLE PRECISION THETGO . THDOT 1 .THD0T2 . DMIn.ETIME 
COMMON/CONST 1/THETGO 110). THD0T1 » THD0T2 .DMIN. ETIME. I YBEG 
C0MM0N/C0NST2/GM . AE » R . RSQ . RQ . THETG . TC 1 
C0MM0N/C0FIT/CS( 5.9) . DA Yl .CENTER .OSTART 

102 FORMAT (1H1.42X33HRUNGE-KUTTA NUMERICAL INTEGRATION/1HO . 48X20HINI 
1AL INPUT VALUES/ 1H020X1HX36X1HY36X1HZ ) 

103 FORMAT (3 (10X. 024,16) ) 

104 FORMAT C5X5HTIME=D24. 16. 2X8HDELTA T=D24. 16. 14X6H0RDER=I2) 

105 FORMAT (lh0.42X33HRUNGE-KUTTA NUMERICAL INTEGRATION/IHO . 45X29HPREDJ 
1CTED EXTRAPOLATED VALUES/1H020X1HX36X1HY36X1HZ) 

106 FORMAT (1H0.42X33HRUNGE-KUTTA NUMERICAL INTEGRATION/lHO .45X29HC0RP E 
1CTED EXTRAPOLATED VALUES/ 1H02QX1HX36X1HY36X1HZ) 

110 FORMAT(1H032HINITIAL COWELL ORDER CHANGED TO 13) 

IF (SW ( 12) ) GO TO 40 
WRITE(3» 102) 

1 WRITE (3. 103) (X(l.I) .1=1.3) . (XU(l.I) . 1=1 » 3) . (XXDO(I) . 1=1.3) 

28 TOUTl=CDEL/TC 
TOUT=T/TC 

WRITE (3. 104) TOUT. T0UT1 .ORDER 
SW ( 1 ) =. TRUE. 

SW(2)=.TRUE. 
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UME=T 

24 KK-K 
2b K.K1 = N 

27 DO 10 J=KK»KK1 
Fj= J 

CSTEPT=TIME+FJ*CDEL 

IF(CSTEPT*TC1+0START .LE. DAY!) GO TO 3 
CALL EPHUAN 
3 K-U+l 
CALL RK 

IF (SW l 12) ) GO TO 10 
IF ( ISWT (8) ) GO TO 11 
ISWT < <3 ) — # TRUE • 

CALL OUTPUT 
ISWT(B)=. FALSE. 

GO TO 10 
11 CALL OUTPUT 
10 CONTINUE 

25 L)0 5 J=1»KK1 
CALL CKDIFF(J) 

5 CONTINUE 

IF(M00E.EQ.4) GO TO 31 

7 SW(3)=.TKUE. 

8 SW(4)=. FALSE. 

SW(b)=. FALSE. 

SW(7)=. FALSE. 

SWl9)=. FALSE. 

CALL TEST 

IF (Sw 19) ) GO TO 31 

30 IF(Sw<4) ) GO TO 29 

1F(SW(5) ) GO TO 16 

IF (SW(7) ) GO TO 2 
GO TO 6 

2 Sw(2)=. FALSE. 

GO TO 24 

16 SW(2)=. FALSE. 

GO TO 8 

29 WHITE ( 3 * 1 1 1 ) 

111 F0KMAT(1H068HINITIAL TABLE FAILED TOLERANCES — PROCEDURE RESTART U)I 
1TH NEW STEPSIZE) 

T=TIME 

K=1 

CALL FRCS 
GO TO 28 

6 IF (SW (2) ) GO TO 31 
wRITEl3#110) ORDER 

31 SW(3)=. FALSE. 

RETURN 

40 IF ( .NOT • ISWT(IS)) GO TO 41 
V»RITE(3» 105) 

GO TO 1 

41 IF (SW ( 14) ) WRITE! 3# 105) 

IF ( .NOT .SW ( 14 ) ) WRITE ( 3 ► 106 ) 

GO TO 1 
ENO 

8 ELT TABLES# 1»671 013# 34119 
0 EOF 8 

SUBROUTINE TABLES 

DOUBLE PRECISION X# XU# XDD » TC # CDEL # XXDD » TOUT # TIME ► TOUT1 » FJ # TT » 
1T»DIFF » TOL1 # T0L2# CSTEPT »FX»FXD»TI»FI1#F C CT #STEPD» TEMP3 
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2»TEMP4,0LDT»TM»FXDD>TREQ 
DOUBLE PRECISION RSAVE » TP ► TPP » TMP 
LOGICAL ISWT > SW 
INTEGER ORDER 

COMMON/TIMING/STEPD (90*3) r TEMP3» TEMP4 >OLDT *TM»XM»YM»ZM»XS» 

* I ENT t ITERS* INCOWL* ICH 

COMMON/RKT/CSTEPT 

COMMON/LIMITS/TOL1 r T0L2 » TC * ISCTf ISWT ( 30 ) » SW ( 20 ) * ORDER » H > L2 . MODE 
COMMON/WORKER/X (200 » 3) ,XD(200»3) , XDD ( 200 r 3 ) t XXDD ( 3) »DIFF<60»3> » 

1 CDEL»TT »T *TREQ»TOUT »K»N 

COMMON/INTERP/FX(60»3) ,FXD(60»3) »FXDD(60*3) »TI »K1»M»M1 
COMMON/TIMES/TABr COWL» FORCS 
COMMON/BRIEF/RSAVE ( 10 ) 

EQUIVALENCE (RSAVE(4) »TP) 

120 FORMAT ( 1H0 * 10X » 22HSTEP-SI2E DOUBLING TO D19.12.6H AT T=.D19,12. 

1 2Xr6H0RDER=.I2) 

122 FORMAT (1H0>10X»22HSTEP-SIZE HALVING TO D19.12r6H AT T=»D19.12. 

1 2Xr6H0RDER=, 12) 

126 FORMAT (lH0»40Xr 33HP0I NTS COMPUTED USING RUNGE-KUTTA) 

CALL CLOCK (XE) 

TOUT=T/TC 
TOUT l=COEL/TC 

C IF OPTIMIZATION OF STEP USE INTERPOLATION 
IF ( ISWT ( 16) ) GO TO 40 
C IF HALVING STEP USE INTERPOLATION 
IF ( .NOT. SW (6) ) GO TO 26 
C IF INSUFFICIENT POINTS USE RUNGE KUTTA 
IF CISCT .LT. 2*(N+1>) GO TO 40 
WRITE ( 3 r 120 ) TOUTlr TOUT » ORDER 
DO 70 1 = 1 » 3 

XDD(K*I)=XDD(K»1)*4.D0 

DO 70 U = 1 t N 

N1=K-J 

C DOUBLE STEP BY SELECTING EVERY OTHER POINT 

N2=K-2*U 
X(Nlf I)=X(N2»I) 

XD(N1»I)=XDCN2»I) 

70 XDD(N1» I)=XDD(N2» I ) * 4. DO 
GO TO 25 

C M ORDER OF HERMITL INTERPOLATION -DEPENDENT ON COWELL ORDER AND STEP 
26 M=( ORDER + 1 ) /2 
M=MAX0(5fM) 

IF ( T EMP3 .GE. 0.5D0) M=MAX0(8rM) 

M=M+2 

IF ( ISCT .LT. M) GO TO 40 
KK=K 

WRITE (3. 122) TOUT 1 r TOUT t ORDER 

K1=K 

M2=M-2 

TEMP3=(-1.D0)*TEMP3 

DO 37 U=2»M2 

FU=U 

C HALVE STEP BY INTERPOLATING FOR MIDPOINTS 
FI1=2.D0*(FU-1.U0) + 1 • DO 
Ti = T-FI 1*CDEL 
M1=M-U 

call hemint 

37 CONTINUE 

TEMP3=(-1.D0)*TEMP3 
M2= M-2 
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DO 30 J1=1»M 
J2=K-J1+1 

XDD( J2» 1)=XDD( J2> 1)/4.DU 
XDD(J2»2)=XDD(J2f2)/4.D0 
30 XDD< J2r3)=XDD< J2.3)/4.D0 
DO 35 J=2»M2 

C INDEX TO MOVE AN ARRAY POINT BACKWARD TO A NEW POSITION 
N1=KK-(M-J) 

C INDEX TO INSERT INTERPOLATED POINT INTO INTERMEDIATE SLOTS IN ARRAY 
N2=KK- 2MM-J) + I 
K=N2 + 1 
DO 32 1=1 * 3 
X(N2»I)=X(N1»I) 

XD (N2* I ) =XD ( Nl » I ) 

XDD (N2» I )=XDD (N1 > I ) 

X ( K » I ) =FX ( J» I ) 

32 XD (K* I ) =FXD( J» I ) 

CALL FRCS 
35 CONTINUE 

C REJECT LAST COMPUTED POINT 
K=KK-1 
T=T-T£MP3 
GO TO 25 

40 TPP=T-2.D0*TEMP3-DBLE(N+1)*CDEL 

C determine if there are enough points at the last STEPSIZE TO PRODUCE 

C N<1 POINTS AT THE NEW STEPSIZE 
IF (TPP.LT.TP) GO TO 50 
ASSIGN 60 TO L 
IF (SW(6) ) GO TO 100 

C IF DECREASING STEP - REJECT LAST COMPUTED POINT 
K=K-1 
T=T-TEMP3 
TPP=TPP-TEMP3 
TT=TEMP3**2 
CALL FRCS 
TT=C0EL**2 
TOUT=T/TC 
GO TO 100 
60 M1=0 

C HERMITE INTERPOLATION ORDER 
M=MAX0(4»ORDER-l) 

IF(M.GT.ISCT) GO TO 50 

C UPPER INOEX OF ARRAY POINTS TO BE INTERPOLATED CK TO K1PI 
K1P=K-ISCT 

C ADJUST TIME SO THAT INTERPOLATION WILL NOT OCCUR NEAR THE END POINT 

C OF THE ARRAY. THESE ARE NOT AS ACCURATELY INTERPOLATED. 

TIME=T-TEMP3 

K1 = K 

C TIME OF NEXT TO LAST POINT OF PARTIAL ARRAY SATISFYING HERMITE ORDER 

C WITHIN WHICH INTERPOLATION MUST OCCUR 

TMP=T-D6LE (M-2) *TEMP3 

C IF TIME OF NEXT TO LAST POINT OF PARTIAL ARRAY OVERLAPS POINTS PRODUCED 

C AT ANOTHER STEPSIZE - SET IT TO WITHIN TPP- TIME OF STEP PRIOR TO CHANGE 

IF(TMP.LT.TPP) TMP=TPP+TEMP3 
KN=N+1 

TEMP3=(-1.D0)*TEMP3 
62 M1=M1+1 
FJ=M1 

C INTERPOLATION TIME 
TI=TIME-FJ*CUEL 
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IF(TI.LT.TMP) 60 TO 65 
CALL HEMINT 

64 GO TO 62 

C UPDATE TIMES FOR NEXT PARTIAL ARRAY 

65 Kl=Kl-(M-3) 

T=T-DBLE(M-3)*DABS(TEMP3) 

TMP=TMP-DBLE ( M-3 ) *D ABS < TEMP3 ) 

IF(KI-M .GE. KIP) GO TO 67 

C IF TIME OF LAST POINT OF PARTIAL ARRAY OVERLAPS STEP CHANGE - ADJU5T TIME 
C AND INDEX TO POINT BEFOR THE CHANGE 
T=T+DBLE (M-K1+K1P-1 ) *DABS ( TEMP3 ) 

TMP=TMP+DBLE ( M-K1+K1P-1 ) *DABS ( TEMP3 ) 

K1=KIP+M-1 

67 M1=M1-1 

IF CTMP.LT.TPP) TMP=TPP+0ABS(TEMP3)-.5D-13 
C IF N<1 POINTS HAVE NOT BEEN PRODUCED CONTINUE INTERPOLATING 
IFCM1.LT. N+l) GO TO 62 
TEMP3=(-1.D0)*TEMP3 

C DROP LAST POINT SINCE INTERPOLATION WAS FROM NEXT TO LAST POINT 
KK=K-1 
T=TIME 
DO 68 1=1. N 
K=KK-I 
DO 66J=1 » 3 

C STORE INTERPOLATED VALUES 
XCK.J )=FX(I.J ) 

66 XD(K»J)=FXD(I»J) 

T=T-CDEL 

CALL FRCS 

68 CONTINUE 
T=TIME 
K=KK 

CALL FRCS 
GO TO 25 

50 IF ( .NOT • ISWT ( 16) ) GO TO 52 
ASSIGN 45 TO L 
IFCSW (6) ) GO TO 100 
ASSIGN 41 TO L 
GO TO 100 

52 IFCSWC6)) GO TO 43 

WRITEC3. 122) TOUT1. TOUT. ORDER 
41 CONTINUE 

C COMPUTE BACKPOINTS USING RUNGE KUTTA 
C 

K=K-1 
T=T-TEMP3 
GO TO 45 
43 FACT=4.D0 

WRITEC3.120) TOUT1.TOUT. ORDER 
XDD(K.1)= XDD(K»1)*FACT 
XDD(K.2)= XDU(K.2)*FACT 
XDD(K»3)= XDD(K»3)*FACT 

45 IF (SW (6) .AND. .NOT . ISWT ( 16) ) GO TO 46 
CALL FRCS 

46 SW(1)=. FALSE. 

WRITE (3. 126) 

TIME=T 

KK=K 

DO 10 U=1.N 
FJ=J 
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cstept=time-fj*cdel 

K=KK-J 
CALL RK 
10 CONTINUE 
K=KK 
T=TIME 

SW ( to) = • TRUE* 

COMPUTE NEW DIFFERENCES AND NEW SUMS 

2S DO 42 J=1»N 
CALL CKDIFF(J) 

42 CONTINUE 
CALL SUMS 
ISCT=N+1 

TP=T-FLOAT(N)*CDEL 
CALL CLOCK (XR) 

TAB = TAB + XR - XE 
RETURN 

C THIS SECTION FORMATS A MESSAGE TO INDICATE THE TYPE OF STEP CHANGE 

100 IF t SW (b) ) GO TO 101 

WRITE ( 3> 121 ) T0UT1 »T0UT » ORDER 

121 FORMAT ( 1H0 r 10X> 21HSTEP SIZE DECREASE TO D19.12* 6H AT T= D19.12 

* 2X t 6H0RDER=r 12) 

GO TO L ( 41 r bO ) 

101 WRITE ( 3 r 123) TOUT1 • TOUT t ORDER 

123 FORMAT ( 1H0 * 10X • 21HSTEP SIZE INCREASE TO D19.12» 6H AT T= D19.12 

* 2Xr bHORDER=» 12) 

GO TO L (45»b0) 

END 

id ELT TEST t l»b71013» 34117 
0 EOF M 

SUBROUTINE TEST 

UOUBLE PRECISION X»XD*XDD»XXUDfDIFF»CDEL>TT >T »TREO 
DOUBLE PRECISION STEPD»TEMP3»TEMP4.0LDT>TM 
DOUBLE PRECISION CSAVE»T0UT»SUM»SUM1>FN1»0PTST 
DOUBLE PRECISION SI ' S2 r PX »PXD»CX» CXD» CTOLf SAVE 
DOUBLE PRECISION T0L1»T0L2>T0L3*T0L4»TC 
LOGICAL ISWT t SW f TEMPO * TEMPI r TEMP2 
INTEGER OROER 

COMMON/WORKER/X(200»3) ► XD C 200 > 3) » XDD (200 • 3) • XXDD ( 3) t DIFF (60 > 3) * 

1 CDEL»TT »T *TREQ»TOUT »K»N 

COMMON/LIMI TS/TOL1 • T0L2» TC» ISCT »ISWT(30)*SW(20) » ORDER » LI »L2» MODE, 
COMMON/TIMING/STEPD (90 » 3) • TEMP3» TEMP4 * OLDT • TM» XM» YM> ZM» XS» 

* I ENT t ITERS » INCOWL ► I CH 
COMMON/COFFS/CSA VE ( 2 ) 

COMMON/OOEL/NL1 » NL2 
COMMON/OPTIM/TOL3»TOL4 

COMMON/COWS/SI (3) »S2(3) »PX(63) rPXD(b3) rCX(63) *CXD(63) >CTOL>SAVE( ItO 
1 » 3) * ITER 

104 FORMAT ( 1H0G4HEST IMATE OF TRUNCATION ERROR INDICATES A COWELL ORDER 
1 CHANGE TO »I3> 9H TIME= » D19.12) 

107 FORMAT (lH0»b4H TRUNCATION ERROR WILL ALLOW AN OPTIMIZATION OF ORD 
*ER CHANGE TO »I3»9H TIME=*D19. 12) 

315 FORMAT ( 1H0 » 47X» 31H0RDER FAILED LI - INCREASE STEP) 

355 FORMAT (1HO#47X»31HORDER FAILED L2 - DECREASE STEP) 

SW(1B)=. FALSE. 

C DIFFERENCE TO BE TESTED 
N1=N-1 
FN1=N1 
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TOUT=T/TC 

IF (SW(3) ) GO TO 5 

C RESTORE POSITION COEFFICIENTS - THE ORDER MAY BE CHANGED 
PX(N+3)=CSAVE(1) 

CX(N+3)=CSAVE(2) 

5 SUM=0.D0 
DO 1 1=1* 3 

1 SUM=SUM + DABS(DIFF(N1. I) ) 

C IF 1ST STEP OPTIMIZATION HAS NOT OCCURRED — DO NOT ALLOW ORDER TO VARY. 
IF ( ( (ISWT(16) ) .AND. SW(20))) GO TO 2 
AFTER MODE 1 STEP CHANGE — ADJUST ORDER 

SWL153 WILL THEN BE SET TO FALSE UNTIL ANOTHER STEP CHANGE OCCURS 
IF ORDER OPTIMIZATION — ORDER OPTIMIZATION CAN OCCUR IF LOCAL ERROR 
SATISFIES TOLLS. 

IF(M0DE.EQ.1.AND.SW(15) ) GO TO 6 
IF ( .NOT • ISWT ( 17) ) GO TO 2 

6 IF (TOL1-SUM) 2,2*7 

7 IF (T0L2-SUM) 6*8*2 

8 IF (Sw ( 3) ) GO TO 2 
SW(15>=. FALSE. 

DO 101 J=1*N1 
N2=N-J 
SUM1=0.D0 
DO 102 1=1*3 

102 SUM1=SUM1+DABS (DIFF (N2* I ) ) 

NN=N2+2 

IF (SUM1.GT • T0L4) GO TO 103 
101 CONTINUE 

SMALLEST ORDER — TO SATISFY THE SPECIFIED TOLERANCE —SIGMA 

103 L0RD=NN+2 
N5=L0RD-3 

IF (LORD.LT , NLl.OR.LORU.GT . NL2) GO TO 2 
IF (LORD .GE. ORDER) GO TO 2 

IF ( (DABS (DIFF (N5, 1 ) ) +DABS (OIFF (N5, 2) ) +DABS (DIFF (N5* 3) ) ) .GT. TOL~ } 

* GO TO 2 

COMPUTE APPROXIMATE LOCAL ERROR FOR NEW ORDER 
SUM=0 .DO 
DO 4 1=1*3 

4 SUM=SUM+DABS(DIFF(NS*I) ) 

ORDER=LORD 
N=0RDER-2 
N1=N-1 
FNl=Nl 

IF (MODE .NE. 1) GO TO 9 
MODE 1 - CHECK LI LESS THAN ORDER LESS THAN L2 
IF ( ORDER, GE.Ll) GO TO 18 
WRITE ( 3 * 107) LORO*TOUT 
WRITE(3*315) 

GO TO 10 

18 IF ( ORDER. LE.L2) GO TO 9 
WRIT E(3*107) LORD, TOUT 
WRITE(3*355) 

GO TO SO 

9 WRITE(3* 107) LORU.TOUT 


2 CONTINUE 

IF (SUM .GE. T0L2 ) GO TO 3 
GO TO (10*20*30) ,MODE 
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C ENTRY MODE 1 - DOUBLE STEP— INCREASE ORDER 

10 IF (ORDER .GT. LI) GO TO 30 
ORDER= LI + (L2-LD/2 
N= ORDER - 2 

C ENTRY MODE 2 - DOUBLE STEP 

20 SW(6)=.TRUE. 

OPTST= (T0L3/SUM*CDEL**FN1 ) ** ( I .D0/FN1 ) 

TEMP0=TEMP1 

TEMP1=TEMP2 

TEMP2=SW(6) 

ISCT1=ISCT2 

ISCT2=ISCT3 

ISCT3=ISCT 

IF ( ICH .LT. 3) GO TO 24 
C TEST FOR INCREASE-DECREASE LOOP. 

IF ( (ISCT1+ISCT2+ISCT3) .GT. 10) GO TO 24 

IF ( (TEMP2. AND. .NOT. TEMPI. AND. TEMPO) .OR. ( . NOT. TEMP2. AND. TEMPI. AND. 

* . NOT. TEMPO) ) GO TO 110 

24 TEMP3=TEMP4 

IF ( ICH.LT . 3) ICH=ICH+1 
CDEL = CDEL * 2.0D0 
TT=CDEL**2 
TEMP4=CDEL 
IF (SW(3) ) GO TO 21 
IF (.NOT. ISWT ( 16) ) GO TO 25 
C STEP OPTIMIZATION — USE COMPUTED STEP SIZE 
CDEL=OPTST 
TEMP4=CDEL 
NSAVE = N 

C SET INDEX TO PRODUCE SUFFICIENT POINTS SO THAT ORDER CAN INCREASE IN MODE 1 
N=L2-2 

IF (.NOT. SW ( 20 ) ) GO TO 25 
SW(20)=. FALSE. 

N=NSAVE 

25 TT= CDEL**2 
CALL CLOCK (YM) 

zm=ym-xm 

C COMPUTE REQUIRED BACKPOINTS 
CALL TABLEB 
SW(1S)=.TRUE. 

IF ( ISWT ( 16) ) N=NSAVE 
CALL CLOCK (XM) 

GO TO 40 

C ENTRY MODE 3 - DECREASE ORDER 

30 ORDER=ORDER-l 
SW(16)=.TRUE. 

IF (ORDER. LT.NL1) GO TO 130 
N=OROER - 2 
IF (SW(3) ) GO TO 22 
WRITE (3.104) ORDER .TOUT 
60 TO 100 
C 

3 IF (SUM .LE. TOLl) GO TO 100 
GO TO (50.60.70) .MODE 

C ENTRY MODE 1 - HALVING STEP — DECREASE ORDER 

50 IF (ORDER .LT. L2) GO TO 70 
0RDER=L1+ (L2-L1 ) /2 
N=0RDER-2 

C ENTRY MODE 2 - HALVE STEP 

60 SW(6)=. FALSE. 
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OPTST= ( T0L3/SUM*CDEL**FN1 ) ** < 1 .D0/FN1 ) 

TEMP0=TEMP1 

TEMP1=TEMP2 

TEMP2=SW(6> 

ISCT1=ISCT2 

ISCT2=ISCT3 

ISCT3=ISCT 

IF ( ICH .LT. 3) GO TO 61 

IF( (ISCT1+ISCT2+ISCT3) ,GT. 10) GO TO 61 

IF ( ( TEMP2 . AND . . NOT . TEMPI . AND. TEMPO > . OR . ( . NOT . TEMP2 . AND . TEMPI . AND . 
* .NOT. TEMPO)) GO TO 110 
61 TEMP3=TEMP4 

IF ( ICH.LT • 3) ICH=ICH+1 

CDEL= CDEL/2.0D0 

TT=CDEL**2 

TEMP4=CDEL 

IF (SW ( 3) ) GO TO 21 

IF( .NOT . ISWT ( 16) ) GO TO 65 

CDEL=OPTST 

TEMP4=CDEL 

NSAVE=N 

N=L2-2 ' 

65 TT= CDEL**2 
CALL CLOCK (YM) 

ZM=YM-XM 

66 CALL TABLED 
Sh (15)=. TRUE. 

67 IF ( ISWT ( 16) ) N=NSAVE 
CALL CLOCK (XM) 

K-K + l 

SW (8) = .TRUE. 

GO TO 40 

C ENTRY MODE 3 - INCREASE ORDER 

70 ORDER=0RDER+l 
SW(18)=.TRUE. 

IF ( ORDER. GT <NL2) GO TO 130 

N= ORDER-2 

IF(SW(3) ) GO TO 23 

IFdSCT.LT. ORDER-1) GO TO 170 

WRITE (3.104) ORDER. TOUT 

K=K-1 

T=T-CD£L 

DO 72 1=1. N1 

DO 72 J=1 . 3 

72 DIFF ( I . J) = SAVE(I.J) 

CALL CKDIFF (N) 

CALL SUMS 
K=K+1 

SW(ti) = .TRUE. 

GO TO 100 
C TERMINATION 

40 TM=T-OLDT 
OLDT=T 
CALL RESUME 
SW(18)=.TRUE. 


100 IF (SW (3) ) RETURN 
CSAV£(l)=PX(N+3) 
CSAVE(2)=CX(N+3) 
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PX(|\I+3)=0.D0 

CX(N+3)=0.D0 

C IF A STEP OR ORDER CHANGE HAS OCCURRED — RETURN 
IF(SW(18) ) RETURN 

c 

C OPTIMIZE STEP FOR FIRST COWELL POINT 
IF ( Sw ( 20 ) ) GO TO 27 
RETURN 

21 SW ( 4 ) =. TRUE • 

I SC T 3=0 
ICH=0 

TEMP2=. FALSE. 

TLMP3=0.D0 

RETURN 

22 IF (ORDER .LE. LI) GO TO 10 
SW(5) =.TRUE. 

RETURN 

23 IF (ORDER .GE. L2) GO TO SO 
SW(7)=.TRUE. 

RETURN 
27 CONTINUE 

C FIRST COWELL POINT -- OPTIMIZE STEP 

OPTST= ( T0L3/SUM*CDEL**FN1 ) ** ( 1 .D0/FN1 ) 

IF (OPTST .LT. CDEL ) GO TO 300 
PX(N+3)=CSAVE(1) 

Cx(N+3)=CSAVE(2) 

IF(M0DE.EQ.2) GO TO 20 
0RDER=Ll+(L2-Ll)/2 
N=OROER-2 
60 TO 20 
300 RETURN 
C ERROR EXIT 

110 Sw(19)=.TRUE. 

M0DE=4 

WRITE ( 3.210) MODE 

210 FORMAT (68H0 POSSIBLE HALVING AND DOUBLING LOOP. MODE WILL BE CHF 
♦NGED TEMP. TO 13) 

IF (.NOT. TEMPI) RETURN 

TEMP3=TEMP4 

CDEL=CDEL/2.D0 

TEMP4=C0EL 

SW(6)=. FALSE. 

GO TO 65 

130 WRITE (3.230) TOUT 

230 FORMAT (1H0. 48H0RDER CYCLE DOES NOT CONVERGE WITHIN SET LIMITS. 

* 1H0 ► 40X . 11HFINAL TIME= 024.12) 

Sw(9)=.TRUE. 

RETURN 

170 0RDER=Ll+(L2-Ll)/2 
N=0RDER-2 
GO TO 60 
END 

Q ELT TNODE. 1 .670830 » 47605 
IS EOF V 

SUBROUTINE TNODE 

DOUBLE PRECISION AT . Z . F J . X . XD . XDD . XXDD , DIFF . CDEL . TT . T . TOL1 . T0L2 . 

1 . TI .PRO .PROl . FX . FXD. A . OUTT . TREQ . TOUT » FXDD 
DOUBLE PRECISION XNODE.XDNODE.TNOD.CC 
DOUBLE PRECISION FND.U. DD . C »D1 . D2 . SUM. SUM1 . ANS . ANSI 
DOUBLE PRECISION ANS2.DD0.D3.SUM2 
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INTEGER ORDER 
LOGICAL ISWT.SW 

DIMENSION D (100*3) *DD(100*3) *DDD(100*3) *01(20*3) *D2(20*3) *D3(20* 
1*ANS(3) *ANS1(3) 

DIMENSION AT (20) * Z( 20 ) , A (20 ) * OUTT (20) 

COMMON/WORKER/X(200*3) *XD(200*3) *XDD(200*3) *XXDD(3) *DIFF(60*3) * 

1 CDEL*TT *T *TREQ*T0UT *K*N 

COMMON/NODE/XNODE (200*3) ► XDNODE (200*3) *TNOD(200) *C(20) *CC(20) * 
IKK ► J1 *ND IFF* NEXT *KEXT * INODE 

COMMON/LIMITS/TOL1 * TOL2* TC ► ISCT* I SWT ( 30) ► S«/ (20 ) ► ORDER ► LI *L2* MODE 

COMMON/INTERP/FX(60*3) *FXD(60*3) *FXDD(60*3) ► TI *K1 *M*M1 

FND=NEXT 

PR0=1.D0 

IN0U£=IN0DE+1 

DO 1 1=1*6 

J=I-1 

L=K-J 

Z( I )=X(L* 3) 

PRO=-Z(I)*PRO 

FJ-J 

AT ( I ) = ( T-FJ*CDEL) 

OUTT ( I ) =AT ( I ) /TC 

1 CONTINUE 
TI=U.D0 

DO 2 J=1 * 8 
PROl=l.D0 
DO 3 1=1*8 
IF(J.EO.I) GO TO 3 
PR01=PR01*(Z(J)-Z(I) ) 

3 CONTINUE 

A ( J) =-PRO/ ( Z ( J ) *PR01 ) 

2 TI=TI+A(U)*AT(J) 

K1=K 

M=8 

Ml=l 

SW(10)=.TRUE. 

CALL HEMINT 
DO A 1=1*3 

XNOD£( INODE* 1)=FX( 1*1) 

4 XUNODE ( INODE* I ) =FXD (1*1) 

TNOD( INODE)=TI 

IF ( ISWT ( 18) ) RETURN 

IF ( INODE. LT.NDIFF) GO TO 10 

DO 11 1=1*3 

DO 11 J=KK * J1 

U2=2+NEXT*(J-1) 

UDD( J* I )= (TNOO( J2)-TN0D( J2-1) ) 

D ( J* I ) = ( XNODE ( J2 * I ) -XNODE ( J2-1 * I ) ) 

11 DD ( J* I ) = (XDNODE ( J2 * I ) -XDNODE ( J2~l * I ) ) 

DO 12 1=1*3 

DO 12 J=1*KEXT 

CALL FKDIFF (D*D1*J*J1*I) 

CALL FKDIFF(DD*D2*J*J1*I) 

CALL FKDIFF(DDD*U3*J*J1*I) 

12 CONTINUE 

IF (SW ( 14) ) GO TO 20 
DO 14 1=1*3 
SUM=0.D0 
SUM1=0.D0 
SUM2=0 • 0U0 
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DO 15 J=1»KEXT 
SUM2=SUM2+C ( J ) *03 ( J » 1 ) 

SUM=SUM+C ( J) *D1 ( J* I ) 

15 SUM1=SUM1+C(J)*D2(J#I) 

J2=(J1-1)*NEXT+1 

ANS ( I ) = ( SUM+U ( J1 » 1 ) ) *FND+XNOUE ( J2 » I ) 

ANSI ( I ) = (SUM1+DD ( Jl » I ) ) ♦FND+XDNODE ( J2» I) 
14 CONTINUE 

ANS2=(SUM2+DDD(J1»1) ) *FND+TNOD ( J2 ) 
INODE=INODE+( NEXT-1) 

T=ANS2 

TNOD ( INOUE ) =ANS2 

DO 16 1 = 1 » 3 

XNODE( INODE r I ) =ANS ( I ) 

XDNODE (INODE * I ) =ANS1 ( I ) 

X(lf I)=ANS(I) 

16 XD(1»I>=ANS1(I) 

X(1»3)=0.0D0 

SW l 12) =• TRUE. 

NDIFF=NDIFF+NEXT 

J1=J1+1 


10 

20 


22 


21 


23 


B ELT 
B EOF 


KK= Jl 

IF ( IS«/T ( 15) ) RETURN 
SW ( 14) = . TRUE • 

CONTINUE 

RETURN 

CONTINUE 

CORRECTING EXTRAPOLATED VALUES 

DO 21 1=1 t 3 

SUM= O.DO 

SUM1=O.DO 

SUM2=O.DO 

DO 22 J=lf KEXT 

SUM2=SUM2+CC ( J ) *D3 ( J » I ) 

SUM=SUM+CC ( J) *D1 (J» I) 
SUM1=SUM1+CC(J)*D2(J#I) 

J2=(Jl-2)*NEXT+l 

X ( 1 » I ) = ( SUM+D ( Jl » I ) ) *FND+XNODE ( J2 1 1 ) 

XD ( 1 » I ) = (SUM1+DD ( Jl 1 1 ) ) *FND+XDNODE ( J2 1 1 ) 
X(l»3)=0.000 

T=(SUM2+DDD( Jl»l) )*FND+TNOD( J2) 

INODE=INODE-l 

J=IN0DE 

TNOD ( J) =T 

DO 23 1=1 f 3 

XNOU£(J»I)=X(l*I) 

XDNODE( J* I ) =XD ( 1 » I ) 

Stt(14)=. FALSE. 

SW(12)=.TRUE. 

RETURN 

END 

TRANS# 1# 670830# 47606 


M 

SUBROUTINE TRANS(TX#TXO) 

DOUBLE PRECISION X# XD # XDD » XXDD#DIFF » CDELr TT » T # TREQ# TOUT 
DOUBLE PRECISION ELEM#TX# TXD#DSlN#DCOSrDSQRT 
DOUBLE PRECISION SING»COSG»SlNH»COSH#SlNI #COSI »SINE#COSE 
DOUBLE PRECISION N#B»P»Q»L#E#DE#C#D#TEP 

CQMMON/WORKER/X (200 #3) »XD(200»3) # XDD (200 #3) »XXDD(3) #DIFF(60#3) # 


1 


CDEL# TT # T# TREQ# TOUT # K # M 


86 


COMMON/EMS/ELEM ( 6 ) *TEP 
DIMENSION TX(3)*TXD(3) 

DIMENSION P ( 3) *Q(3) 

SING=DSIN(ELEM<5)) 

COSG=DCOS ( ELEM ( 5 ) ) 

SINH=DSIN(ELEM(6) ) 

C0SH=DC0S(ELEM(6) ) 

SINI=DSIN(ELEM(3) ) 

COSI=DCOS(ELEM(3) ) 

N=1.0D0/(DSQRT(£LEM(1)**3) ) 

B=ELEM ( 1 ) *DSfclRT ( 1 . ODO-ELEM ( 2 ) **2 ) 

P ( 1 ) =cosg*cosh-sing*sinh*cosi 

P ( 2 ) =SING*COSH*COSI+COSG*SINH 
P(3)=SING*SINI 

Q ( 1 ) =- ( SING*COSH+COSG*SlNH*COSI ) 

Q ( 2) =COSG*COSH*COSI-SlNG*SINN 
Q(3)=C0SG*SINI 

2 L=ELEM ( 4 ) +N* ( T-TEP ) 

L=DMOD (L * 6. 2831853071795864D0 ) 

E=L 

3 DE= L-E+ELEM(2)*DSIN(E) 

IF<DABS(DE)-0.2D-14.LE.0.0D0) GO TO 5 
E=E+D£/ ( 1 . ODO-ELEM (2) *DCOS (E+0 . 5D0*DE > ) 

GO TO 3 

5 SINE= DSIN(E) 

COSE=DCOS(E) 

C=C0SE-ELEM(2) 

D=1 . ODO-ELEM ( 2 ) *COSE 
D=N/D 

DO 10 1=1*3 

TX ( I ) =ELEM ( 1 ) *P ( I ) *C+ B*Q(I)*SINE 
10 TXD(I) =D* (B*COSE*Q( I ) -ELEM ( 1 ) *SINE*P( I ) ) 

RETURN 

END 

8 ELT TRANS2* 1 * 670630 * 47607 
8 EOF 8 

SUBROUTINE TRANS2 (X*XD) 

C TO COMPUTE ELEMENTS FROM POSITION AND VELOCITY VECTORS 

DOUBLE PRECISION X*XD*ELEM*RSQ*VSQ*RRD*P.PSQ*DATAN2*DSQRT*ES 

iine*ecose*r*e*y*datan*sinh*cosh*sini*dsin*q*sinu»sinv*cosu.cosv 

2*TEP 

DIMENSION X(3) *XD(3) *P(3) 

C0MM0N/EMS/ELEM16) *TEP 

RSQ=O.DO 

VSQ=O.DO 

RRD=O.DO 

PSGi=O.DO 

P ( 1 ) =X (2) *XD (3) -X ( 3) *XD (2) 

P (2) =X ( 3) *XD ( 1 ) -X ( 1 ) *XD (3) 

P(3)=X<1)*XD(2)-X(2)*XD(1> 

DO 1 1=1*3 
RSQ=RSQ+X ( I ) **2 
VSQ=VSG+XD ( I ) **Z 
RRD=RRD+X(I)*XD(I) 

1 PSQ=PSQ+P ( I ) **2 
Q=DSQRT (PSQ) 

R=DSURT (RSQ) 

ELEM ( 1 ) = R/(2.D0-R*VSQ) 

SINV= Q*RRD/R 
COSV= (PSQ/R)-1.D0 
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ELEM(2)=DSQRT(SINV**2+C0SV**2) 

SINV=SINV/ELEM(2) 

C0SV=C0SV/ELEM(2) 

ESINE=R*SINV/ (ELEM ( 1 ) *DSQRT ( 1 .DO-ELEM (2) **2) ) 
ECOSE=R*COSV/ELEM ( 1 ) +ELEM ( 2 ) 

E=DATAN2(ESINE»EC0SE) 

ELEM ( 4 ) =E-ELEM ( 2 ) *ESlNE 
r=DSQRT(1.0D0-(P(3>/O)**2)/(P(3)/Q) 

ELEM(3)=DATAN(Y) 

SINI=DSIN(ELEM(3) ) 

SINH=P(1)/(Q*SINI) 

C0SH=-P(2)/(G*SINI) 

ELEM ( 6 ) =DATAN2 ( SI NH r COSH ) 

SINU=X(3)/(R*SINI) 

COSU= ( X ( 1 ) *COSH+X ( 2 ) *SI NH ) /R 

ELEM ( 5 ) =DATAN2 ( SINU » COSU ) -DATAN2 ( SI NV * COSV ) 

IF(ELEM(4> .LT.O.ODO) ELEM(4)=ELEM(4)+6.2831853071795864D0 
IF (ELEM(5) .LT.O.ODO) ELEM ( 5) =ELEM ( 5) +6. 2831853071795864D0 
IF (ELEM(6) .LT.O.ODO) ELEM(6)=ELEM(6)+6.2831853071795864D0 
RETURN 
END 

0 ELT YMDAY>lr 670830 1 47607 
8 EOF 8 

DOUBLE PRECISION FUNCTION YMDAY ( I YMD» iHMr SEC ) 

DOUBLE PRECISION SEC 
IY=(IYMD/10000) *10000+101 
1HMS=IHM *100 

CALL DIFF ( I Y * 0 r I YMD »IHMS»IDrIS) 

YMDAY=86400*(ID+1)+IS 
YMDAY= ( YMDAY+SEC ) /8. 64D4 
RETURN 
END 


8 ELT EPHQANr 1 r 670913 r 46904 
8 EOF 8 

SUBROUTINE EPHQAN 

DOUBLE PRECISION DUMMY r C r DAY1 r CENTER* D. TIME ( 1 1 )» VAR ( 11 ) 
DOUbLE PRECISION FI 
DIMENSION AO (9. 11) 

DOUBLE PRECISION DSTART 

C0MM0N/C0NST1/DUMMY(14) rlYBEU 

COMMON/COF IT/C ( 5 . 9 )* DA Y1 . CENTER r DSTART 

IYl=IYBEG-59 

CENTERED AY 1+«J . 5D0 

DO 5 I=l»ll 

F1=I-1 

D=DAY1+.5D0*F1 
TIME(I)=D-CEnTER 
5 CALL EPHEM (IYlrDrAO(lrl) ) 

DO 10 I— 1 > 9 
DO lb J— 1 » 11 
15 VAR(J)=AO(I»J) 

10 CALL COEFF(VAR* TIME »11»4»C(1»I)) 

DAY1=DAY1+5.D0 

RETURN 

END 
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APPENDIX A 


RUNGE KUTTA INTEGRATION FORMULAS 
K 0 = hf(X 0 ,Y 0 ) (i) 

K, = h f(x„ + ^-h, Y„ t^Ko) ( 2 ) 

K 2 = h f (x„ +yh, Y„ + ^ <K 0 * 3K,)) (3) 

^3 — ^ f + "X^’ Y 0 + (K 0 + 3K 2 )^ (4) 

K 4 = h f (x 0 + y h > Y o + y( K o + 3K 3 )) (5) 

K s = hf (x 0 + yh, Y 0 + (13K 0 - 27K 2 + 42K 3 + 8 K 4 )j ( 6 ) 

K 6 = h f + yh, Y 0 + 43 ^- (389K 0 - 54 K 2 + 966 K 3 - 824 K 4 + 243K S )) (7) 

K 7 = hfj^+h, Y 0 +^ (-231K 0 + 81K 2 - 1164K 3 

+ 656 K 4 - 122K 5 + 800K 6 )j ( 8 ) 

K 8 = h f (x 0 +|h, Y 0 + 2 §g (" 127K o + 18K 2 - 678K 3 

+ 456 K 4 - 9 K g + 576K 6 + 4K 7 )] (9) 
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K g = h f |x 0 + h, Y 0 + (148lK Q - 81K 2 + 7104K 3 - 3376K 4 

+ 72K S - 5040 K 6 - 60K 7 + 720K g )j (10) 

Y (X Q + h) = Y(X 0 ) + (41K 0 + 27K 3 + 272K 4 + 27K S 

+ 216 K 6 + 216 K g + 41 K 9 ) 


Reference: Shanks, E. (1966): "Solutions of Differential Equations by Evaluations of Functions," 
Math, of Compnt., Vol. 20, pp 21-38. 
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APPENDIX B 


FORCE MODEL 


The equations of motion used in the program are of the form 


X 


+ X < + x 0 + ^ + X, 


where 

(A) X e - the earth’s gravitational effects 

(B) X f , X q - lunar and solar gravitational effects 

(C) Xjj - earth's atmospheric drag effects 

(D) X, - solar radiation effects. 

The computation of the required earth, lunar and solar ephmeris quantities re- 
quired by these perturbative accelerations is described in Section (E). 

In what follows, the formulation used to compute X(t), given X(t), X(t) and t 
is detailed. We remark that the unit of length = 6378.166 km and unit of time = 
806.81242 secs. 

(A) EARTH’S GRAVITATIONAL PERTURBATION 
The earth's gravitational perturbation is given by 




/xX d/x 

r 3 BX 


where U, the disturbing potential, expressed in the spherical coordinates of the 
satellite r, \p, is given by 


Reference for section (A) and (E): WRDC Final Technical Report, "Unified Geodetic Parameter 
Program (GEOPS)”, December 1964 to March 1966, Vol. 1, Mathematical Analysis by Kahler and 
Wells, pp 24-28, appendices A and B. 
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u 


ii 

r 


L LW 


n = 2 m = 0 


(C cos mk + S sin mX.) P® (sin i p) 

v nm nm / n x ~ ' 


where 

b = mean radius of earth (unit of length) 
r = geocentric radius of satellite 

A. = longitude of satellite measured eastward, given by A. = a - 9 g , 

where a, the right ascension is given by a = tan -1 (y/x), and 6 g is 
the right ascension of Greenwich at time t , referenced to true 
equinox (see Section E) . 

Sin tp = sine of the geocentric latitude, given by Sin \p = z/r 

C , S = the unnormalized harmonic coefficients (stored as data arrays 
with a maximum index of 20) . 

P“ = the m th derivative of the n th Legendre polynomial 

fj , = product of gravitational constant and mass of the earth. 

The perturbative force, in rectangular coordinates is given by 

9U 3U3r_ jHjBA^ 3U 30 
3X ~ 3X + BX 30 3X 

where 


3U 
3 r 


J±_ 

v 


N n 

V ' V — l 


b n (n + 1) 

n=l m*0 


ICC m + SS ) 

l nm nm J 


-n + 1 


pm 

n 


(Sin <//) 


92 



P n m (Sin 4 ,) 


BU 

BX 


£ 

r 


N n 



n=0 m = 0 


{sc nm - cs \ 

l nm nmj 


N n 


BU _ (± 
Bi// r 


E E {p "‘ 


n = 0 m = 0 


(Sin \p) - m tan \p P“ (Sin <//)} 


(CC + SS ) 

V. nm nmj 


where 


CC 


SS 


sc 


cs 


C Cos mX 

nm 

S Sin mX 

n m 

S Cos mX 

nm 

C Sin m\ 

nm 


and are computed using the relations 

Cos mX = 2 Cos X Cos (m - 1)X - Cos (m - 1)X 

Sin mX = 2 Cos X Sin (m - 1)X - Sin(m - 1)X 
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and where the P m 

n 


(Sin i p) are given recursively by 
p n° = [( 2n - 1) Sin 'P P n— i - (n - 1) p n °_ 2 ] /n (m = 0) 

P n = ( 2n ~ 1) Cos \p P^Jj 1 (m = n) 

p n = P n-2 + ( 2n - 1) Cos ^ P™;* (m J 0, m ¥ n) 

Finally, the position partials of r, X, \p are given by 

B r _ /_x _y _z_\ 

BX ~ \ r ’ r ’ r ) 

BX _ / BX BX BX \ 

BX \ Bx ’ By ’ Bz ) 

where 

BX _ ~ y BX _ x BX 

Bx x 2 + y 2’ By x 2 + y 2’ Bz " ° 

30 _ / Bi/» B \p 3 0 \ 

BX \ Bx ’ By ’ Bz / 

where 

B0 _ - xz Bi/j _ ~ yz 30 _ Tx 2 + y 2 

Bx 2 V 2 ! 2 ' By 2 2 i 2 ' Bz 2 

r* Tx^ + y z J i"* rx 2 + y* r 

(B) LUNAR AND SOLAR GRAVITATIONAL EFFECTS 
The moon's gravitational perturbation is given by 
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where 



X = satellite position vector 

\ = moon's position vector expressed in the geocentric system (Section E) 

= the lunar mass expressed in earth masses 
r ( g = satellite-! moon radius vector 
^ = moon radius vector, and 

P = product of gravitational constant and mass of the earth. 

The solar gravitational effect is computed precisely as above, with X <t M , 

*.» r c replaced by X 0 , M # , r 0s , r 0 . 

(C) ATMOSPHERIC DRAG EFFECTS 

The drag effect on the accelerations is given by 

„ Cd a 

*DRAG ~ 2M + ^l] [* + <°2 t ] /°(^) l^r^r 

where 

Cp = atmospheric drag constant 

A = cross-sectional area of satellite in CM 2 

M = mass of satellite in gins. 

P x ,p 2 - atmospheric density correction parameters 

/O(h) = atmospheric density at satellite height h 

V f = relative velocity vector of the satellite, given by (x - w e y, y - w x, z), 
where w e is the angular velocity of the earth in rad./c.u.t. 

9 = the angle denoting the lag between the atmospheric bulge angle and 
local noon. 
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The atmospheric density p(h) is given by 


P(h) = P N ( h) + [p 0 ( h) - p N (h)] Cos'" (0/2) 

where p N (h), p D (h) are night and day density values computed from a set of 
tabular values p N (hj), p D (hj) as follows: 

A table look-up is performed to obtain the index i satisfying 

h. < h < h. +1 , 

1 1 + 1 ’ 


from which p N (h) or Pj) (h) are computed by 


Pj ( h ) 


Pi ( h i ) 


xp 



- h 
H. 


»j = h i -hi M /ln (h^] 

where j = 1 for night density, and j = 2 for day density. 

Also, Cos 0 is given by 

Cos 0 = X* • X* ulge 
where * denotes unit vector and 

X* = ( X* Cos 6 - Y * Sin 6 , Y * Cos 6 + X* Sin 6,1*) 

where (X*, Y*, Z*) is the unit vector of the sun’s position. 

' © © © 

(D) SOLAR RADIATION EFFECTS 

The effects due to solar radiation are given by 
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where 


crA ( X e~ X) 
M |X 0 - X| 


A = cross-sectional area of satellite in cm 2 

M - mass of satellite in gins 

X q = sun's position vector 

a - solar radiation constant dyne/cm 2 

X = satellite position vector 

The sunlight -shadow determination is given by 

X*.X* > 0 = the satellite is on the solar side of a plane normal to the 
earth-sun line. No shadowing can occur. 

X* . X * < 0 = the satellite will be shadowed only if the distance between 
the satellite and the earth-sun line is less than the earth's 
radius, hence if 

|X x X*| < T 
© 

where 

T = 1 - f (z + Z @ /r 2 - 1 ) 2 , 

then the satellite is in the earth's shadow. Otherwise, it is 
in sunlight. 

Note: * denotes unit vector, and f is the flattening coefficient of the earth. 


(E) EPHEMERIS QUANTITIES 

The right ascension of Greenwich at the time of interest (t) is com- 
puted from 


97 


where 


0% - £g 0 + (t - t 0 ) e x + (t - t 0 ) e 2 + A/i 


6g 0 = right ascension of Greenwich at Jan 0.0 for the year of interest 

6 1 - 0.9856473354 deg/mean solar day 

6 2 = 360.9856473354 deg/mean solar day 
A// = equation of equinoxes 

(t - t 0 ) = time in days from Jan 0.0 for the year of interest. 

The equation of equinoxes A/i is obtained in the process of computing the lunar 
and solar ephemerides. 


The ephemeris quantities are computed at ten equal intervals over a five 
day period and least squares fit to a fourth order polynomial in time about the 
midpoint of the five day period. The positions are determined for intermediate 
points by evaluating the polynomial at the required time. 


To calculate the earth-centered slant range and inertial unit vectors to the 
sun and moon and the equation of the equinoxes at a given time, the time interval 
from epoch (1900 Jan 0.5) is denoted by T when measured in Julian centuries, 
by D = 3.6525T and d = 10000D. 


Angular variables: 


e - mean obliquity of the ecliptic 
L & - geometric mean longitude of the sun, mean equinox of date 
T - mean longitude of solar perigee 

- mean longitude of the moon, measured in the ecliptic from the 
mean equinox of date to the mean ascending node of the lunar 
orbit, and then along the orbit 

T'- the mean longitude of lunar perigee, measured in the ecliptic 
from the mean equinox of date to the mean ascending node of the 
lunar orbit, and then along the orbit 
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n 


V 


F 


D 

Ay 

Ae 


e 
A jj. 


X , Y , Z 
© ©’ © 


WZic 


X. 


<c 


the longitude of the mean ascending node of the lunar orbit on the 
ecliptic; measured from the mean equinox of date. 

L c minus the mean longitude of the lunar perigee, measured in 
the ecliptic from the mean equinox of date to the mean ascending 
node of the lunar orbit, and then along the orbit 

geometric mean longitude of the sun minus the mean longitude of 
perigee of the sun 

L c minus the longitude of the mean ascending node of the lunar 
orbit on the ecliptic; measured from the mean equinox of date 

L c minus the geometric mean longitude of the sun. 

nutation in longitude 

nutation in obliquity 

true obliguity of the ecliptic 

equation of equinoxes or nutation in right ascension 

eccentricity of sun in an earth centered coordinate system 

semi major axis of earth's orbit about the sun (cul) 

radius vector to sun 

apparent longitude of sun 

inertial unit vectors to sun 

semi major axis of lunar orbit (cul) 

radius vector to moon 

eccentricity of lunar orbit 

inertial unit vectors to moon 

apparent longitude of moon 


99 


4> “ apparent latitude of moon 

v - true anomaly of the sun at Jan 0.0 for the year of interest 
6 x - deg/mean solar day. 

e = 23°. 452294 - 0°. 0130125 T - 0.164 x 10" 5 T 2 

L = 279° 41 ' 48" . 04 + .12960276813" x 10 9 T + l" .089T 2 

© 

T = 281° 13' 15" .00 + 6189" .03T + l" .63T 2 

L, = 270°. 434164 + 13° . 1763865268 d - .85 x 10" 4 D 2 

T' = 334°. 329556 + 0° . 1114040803 d - .7739 x 10" 3 D 2 
n = 259°. 183275 - 0° . 0529539222 d + .1557 x 10" 3 D 2 
£ = 4 - T' 

£' = L - T 

0 

f = 4 - o 

D = L - L 

c © 

Ay h .2088 sin 20- (17". 2327 + .01737T) sin Q 

- 1.273 sin (F - D + 0) - .2037 sin (2F + 20) 

+ .1259 sin (£') - .0496 sin (£ ' + 2F - 2D + 20) 

+ .0214 sin (-£ 1 + 2F - 2D + 20) + .0675 sin (£) 

- .0342 sin (2F + 0) - .0261 sin (£ + 2F + 20) 

+ .0114 sin (-£ + 2F + 20) + .0124 sin (2F - 2D + 0) 

- .0149 sin (£ - 2D) 
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Ae = 9". 2106 cos 0 - .0904 cos 20 

+ .552 cos 2(F - D + 0) + .0884 cos 2 (F + 0) 

+ .0183 cos ( 2 F + 0) + .0216 cos (£' + 2F - 2D + 20) 

e = e 0 + Ae 

A/x = Avj/ cos e 

e = .01675104 - .418 x 10' 4 T 

© 

£ © = L © + 2e © sin <*1 (t " t o> ' v ) 

X © = cos ©> 

Y © = sin (* 0 ) cos e 

Z © = sin ( £ ©> sin 6 

= % o - ‘V4 + - T )] 

\ = [ 206265 L + 22640 sin(£) - 4586 sin(£ - 2D) 

+ 2370 sin (2D) + 769 sin (2£) - 668 sin (£ ’ ) 

- 412 sin (2F) - 212 sin (21 - 2D) 

- 206 sin (£ + £' - 2D) + 192 sin(£ + 2D) 

- 165 sin (£ ' - 2D) + 148 sin (£ - £') 

- 125 sin (D) - 109 sin (£+£')- 55 sin (2F - 2D) 

- 45 sin (£ + 2F) + 40 sin(£ - 2F) - 38 sin (£ - 4D) 

+ 36 sin (3£) - 31 sin(2£ - 4D) + 28 sin(£ - £' - 2D) 

- 24 sin (£ ' + 2D) + 19 sin(£ - D) 

+ 18 sin (£ ' + D)] /206265 
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18520 sin (S) {l - .00293 sin (1.403808 

- . 0009242203 T)} - 31 sin(F - £ - 2D) 

- 25 sin (F - 2£ ) - 23 sin (£ 1 + F - 2D) 
+ 21 sin (F - £) - 526 sin(F - 2D) 

+ 44 sin (£ + F - 2D)J /206265 

COS 4 > cos X 

cos <p sin X cos e — sin <£> sin e 
cos cp sin X sin e + sin <p cos e 

( 1 - e^. 2 ) / ( 1 + e ( cos £ ) 
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